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PREFACE 


This  is  Volume  II  of  a three-part  report  which  describes  a computer  program 
for  the  prediction  of  Solid  Propellant  Rocket  Motor  Performance.  The  computer 
program  described  herein  will  be  refered  to  as  the  SPP  program. 

Volume  I of  this  report  describes  the  engineering  analysis  which  was  used 
in  developing  this  computer  program. 

Volume  II  of  this  report  is  a programming  document  of  the  computer  program 
which  was  developed  under  this  contract.  It  includes  a subroutine-by  subroutine 
description  of  all  o:c  the  elements  of  the  SPP  program. 

Volume  III  of  this  report  is  a Program  User's  Manual  which  describes  the 
input  necessary  to  execute  the  SPP  computer  program  and  the  information  required 
to  interpret  the  output.  A sample  case  is  also  included. 
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1.0  Introduction 

This  document,  Volume  II,  provides  a detailed  description  of  all  of  the 
elements  used  in  the  Solid  Propellant  Rocket  Motor  Performance  Prediction  (SPP) 
computer  program.  The  goal  in  developing  the  SPP  code  was  to  develop  a computer 
ized  analysis  technique  that  could  predict  solL’  propellant  rocket  motor  delivered 
specific  impulse  to  within  +2%  and  delivered  thrust  and  total  impulse  to  within 

ts%. 

The  contents  of  this  volunu  are  not  intended  to  supply  the  reader  with  a 
complete  description  of  the  SPP  computer  code.  Instead,  it  is  intended  to  give, 
along  with  the  information  supplied  in  Volumes  I and  III,  the  necessary  information 
about  the  internal  functions  of  the  SPP  code  to  allow  the  educated  user  to  under- 
stand and/or  mod  if  v the  basic  internal  workings  of  this  computer  program. 

The  SPP  code  consists  of  five  basic  computational  modules  plus  a master 
control  module  which  selects  via  user  input  which  of  the  computational  modules 
are  to  be  executed.  The  approach  used  in  developing  this  code  has  been  to 
supply  the  user  with  the  maximum  amount  of  flexibility  consistent  with  ease  of 
use. 

The  five  basic  computational  modules  used  in  the  SPP  code  are  described 
very  briefly  below  it.  Table  1-1.  A more  detailed  description  of  these  modules  is 
presented  in  Volume  T and  to  a lesser  extent  in  Volume  III. 


Module 

Description 

0DE 

Calculates  the  reference  or  Realized  performance 
of  the  solid  propellant  rocket  motor. 

BAL 

Calculates  the  grain  shape  and  internal  ballistics 
parameters  as  a function  of  time. 

0DK 

Calculates  the  nozzle  performance  considering  the 
degradation  in  I due  to  finite  rate  chemical  kinetics 

Sp 

TD2P 

Calculates  the  nozzle  performance  considering  the 
degradation  in  performance  due  to  two  dimensional 
and  two  phase  flow  effects. 

TBL 

Calculates  the  nozzle  performance  considering  the 
degradation  in  performance  due  to  a viscous  turbulent 
boundary  layer. 

1-1 


Volume  I consists  of  the  engineering  analysis  incorporated  in  the  SPP  code 
while  Volume  III  describes  the  input  necessary  to  execute  the  SPP  code.  The 
user's  manual  (Volume  III)  and  this  volume  are  considered  to  be  the  prerequisite 
amount  of  informat.on  needed  to  successfully  understand  and/or  modify  the  SPP 
program. 

Section  2 of  this  volume  describes  the  overall  flow  of  the  SPP  program  along 
with  the  data  communicated  between  the  various  computational  modules.  The  over- 
lay structure  and  labeled  common  blocks  used  by  each  subprogram  are  detailed  in 
section  3 while  the  linkages  and  files  used  by  the  SPP  code  are  described  in  sec- 
tion 4.  A complete  subroutine  by  subroutine  description  of  this  code  is  included 
in  section  5.  Section  5 is  orgainzed  by  overlay  structure  and  includes  a description 
of  the  main  routine  of  each  overlay  followed  by  a description  of  each  subprogram 
in  that  overlay  in  alphabetical  order. 

Liberal  use  of  existing  documentation,  References  1 through  6,  about  the 
'ndividual  elements  of  the  SPP  code  has  been  used  in  this  volume.  The  interes- 
ted reader  would  be  well  advised  to  obtain  these  documents  if  he  does  not  already 
have  access  to  them. 


2,0  SPP  Program  Flow  Charts 

This  section  cf  Volume  II,  the  Computer  Program  Description,  is  intended 
to  give  the  user  of  the  SPP  code  a general  description  of  the  overall  flow  of  the 
computer  program.  Also  described  here  is  how  data  is  transmitted  between  modules 
when  the  many  options  which  are  available  in  the  SPP  code  are  executed. 

Figure  2-1  shows  the  overall  logical  flow  of  the  SPP  code.  The  solid  lines 
in  this  figure  show  the  mandatory  flow  while  the  dashed  lines  indicate  the  optional 
paths  the  computer  program  can  take.  In  the  case  that  more  than  one  optional 
path  can  be  taken,  a diamond  sign,0  , is  used  to  indicate  the  alternate  paths. 

The  boxes  on  the  right  hand  side  of  Figure  2-1  show  the  action  taken  by 
the  SPP  code  in  response  to  Input  directive  cards  and  Iso  to  the  selection  of 
loss  modules  to  be  executed.  The  items  on  the  left  band  side  of  this  figure  show 
the  alternate  logical  execution  paths  which  are  determined  by  user  input. 

Figure  2-2  shows  the  actual  execution  sequence  taken  by  the  SPP  code. 

The  following  notes  are  intended  to  clarify  Figure  2-2. 

Alternate  linkage  data  read  in  under  the  control  of  subroutine  PR0B 
is  read  in  in  sequence  (See  Volume  III,  Section  2.11). 

Linkage  data  written  by  the  various  loss  modules  is  positioned  by  sub- 
routine SETUP  before  being  read  by  subroutine  PICKUP.  This  insures 
that  none  of  the  linkage  data  is  overwritten  unless  multiple  cases  are 
executed  (in  which  case,  the  linkage  data  for  the  first  case  is  destroyed). 


2-1 


2-3 


3,0  Overlay  ani  Common  Block  Description 


Section  3,1  describes  the  overlay  structure  of  the  SPP  while  Section  3.2 
describes  the  common  blocks  used  in  the  SPP  code. 


3,1  Program  Overlay  Structure 

The  SPP  Program  overlay  structure  is  shown  in  Table  3-1  ai\d  corresponds 
to  the  order  of  appearance  of  the  subroutine  write  ups  in  Section  5,  if  read  down- 
ward and  then  to  the  right. 

This  overlay  structure  was  developed  under  the  restrictions  of  the  CDC 
6000  series  loader  (SCOPE  3.3)  which  allows  only  three  levels  of  overlay.  Con- 
siderable savings  in  central  memory  storage  can  be  achieved  with  a more  exten- 
sive overlay  structure  is  used. 
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BLOCK  DATA 
ATMOS 


Table  3-1  Program  Overlay  Structure 
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3.2  Common  Blocks 

The  common  blocks  and  the  subroutines  that  are  referred  to  by,  are  listed 
in  Table  3-2.  Tables  3-3  through  3-14  list  the  variables  in  some  of  the  more  im- 
portant labeled  common  blocks  used  in  the  SPP  code.  These  tables  are  organized 
alphabetically  by  common  block  name. 

Appendix  B of  Reference  1 contains  a dictionary  of  the  variables  appearing 
in  labeled  common  blocks  in  the  0DE  module.  While  there  are  somejninor  vari- 
ances between  Reference  1 and  the  SPP  code,  this  dictionary  should  be  very  use- 
ful to  anyone  attempting  to  modify  the  ODE  module.  Unfortunately,  similiar 
dictionaries  are  not  available  for  the  other  modules. 


Table  3-2  Common  Block  References 


COMMON 

BLOCK 

REFERRED  TO 
BY  SUBROUTINE 

ADD 

BLKOAT 

CO'JTR* 

figi 

FIG* 

LOOPS 

PRISM? 

UNION 

amachr 

FAM 

TBLEDG 

ANMft 

PRINT 

BIKDAT 

DRIVER 

ODES 

UOKINP 

OUTPUT 

OUT  1 

PACK 

REACT 

RKTOUT 

SEARCH 

ARPRNT 

NZGEQP 

ODES 

ookinp 

PICKUP 
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COMMON 

BLOCK 


REFERRED  TO 
BY  SUBROUTINE 


BAUOH 

BALL 

CARD 

CDELMX 

cniF 

CIDAU 

CINT 

CK 

CQOATA 


SUMRY 

AVEBAL 

NZCF.Oh 

PICKUP 

SURRY 

BALLS 

SLKDAT 

LOOPS 

MAINJ 

PVS 

READIN 

READIN 

STORE 

DRIVER 

ODES 

SAVE 

START 

zetait 

BARCON 

BARPRO 

BARSET 

getpt 

QUITS 

READSI 

SEVAL 

START 

IAUX 
INT 
LESk 
M A I N 1 D 

ERROR 

READIN 

BARCON 

BARPRO 

barset 

EDI  TCP 

GETPT 

QUITS 

READSI 

START 

zetait 


COMMON 

BLOCK 

REFERRED  TO 
BY  SUBROUTINE 

CODE 

(U.KDAT 

ERROR 

READIN 

STORE 

COOTOK 

ELU 

C0EF5X 
COF I IF 

COHCAS 


COMERS 

COMXP 

CQMy 


MA1N1D 

DDK 

ODKIhP 

CPHS 

SEARCH 

BARPRO 

PIU 

QUITS 

CONVRT 

OERIV 

DRIVER 

EF 

flu 

gtf 

IAUX 

INT 

maimd 

OPKINM 

OUTPUT 

PACK 

PRES 

PEAXIN 

SELECT 

STf 

T T APE 

DRIVER 

ODES 

CONvRr 

DERIV 

DRIVER 

EF 

*LU 

GTF 

IAUX 

InT 

MAIN1D 
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COMMON  REFERRED  TO 

BLOCK  BY  SUBROUTINE 


CONI 

COM2 

COM  J 
COM4 
C0M5 

CPEfP 

CPRNT 

C3EVAL 

C3PRXC 


ODKINM 

OUTPUT 

PACK 

PRE3 

»EAXIN 

STP 

BLKDAT 

TAB 

TABIN 

BLKDAT 

MAINJ 

READJN 


BLKDAT 

MAINS 

READIN 

BLKDAT 

TAB 

TABIN 

BLKDAT 

DATAIN 

READIN 

TAB 

TABIN 

CPH3 

ROCKET 

DRIVER 

INT 

MAIN1D 

OOKiNP 

PACK 

BARPRO 

BARSET 

GETPT 

QUITS 

READSI 

SF.VAL 

START 

CONVRT 

OERIV 

If 

PACK 
RE  AX  IN 
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COMMON 

REFERRED  TO 

BLOCK 

BY  SUBROUTINE 

SELECT 

STOICC 

cutil  cn^vRT 

CPH3 

DRIVER 

FLU 

IAUX 

HAIN1D 

MUK 

odes 

OOK 

UOKINP 

OUTPUT 

OUT  1 

PACK 

PICKUP 

PRES 

REACT 

REAXIN 

PKTUUT 

ROCKET 

SEARCH 

select 

SPLN 
STF 
SUHRY 
TT  APE 

c M all  o»ivm 

FL'J 

« a : m i n 

OHK INP 

PACK 

PRES 


OOLOUP 

DEPIV 

FF 

I AUX 

INT 

MAINJO 

flDKlNP  t 

PRNTCK 

DOUBLE 

EDLHPM 

GAUSS 

MATRIX 

ODES 

OUT1 

E3GEX 

RARPPO 
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COMMON 

BLOCK 


REFERRED  TO 
BY  SUBROUTINE 


EDGCV 

m 

ERR 

EXTR^P 

FIG 

FnALBT 

MHTCOM 

1001 


REAOfI 

BARPfiQ 

GETPT 

REAOSJ 

PROBI.H 

SUHRV 

AD  JK 

AGP 

CNTRt 

ERROR 

FCA'.C 

JArtES 

KPBPT 

N?HAIN 

PART II 

TD2P 

WALL 

«LPT 

DRIVER 

FIND 

PICKUP 

SURRY 


BLKDAT 

CONTRi 

CONTR2 

FIG  1 

FIG2 

PRISMl 

PRI3M2 

READIN 

STORE 

DERI  V 
flu 

IAUX 

HAIN1D 

EQLBRM 

ODES 

AGP 

EISP 

ISP1D 


COMMON 

BLOCK 

REFERRED  TO 
BY  SUBROUTINE 

1N0X 

CPHS 

DRIVER 

EOLBRH 

FROZEN 

GAUSS 

hcalc 

MATRIX 

ODES 

OUTt 

REACT 

RKTOUT 

ROCKET 

SAVE 

SEARCH 

iNt'XX 

CPHS 

EQLBRM 

\ 

FROZEN 

GAUSS 

HCALC 

MATRIX 

OOFS 

REACT 

RKTOUT 

ROCKET 

SAVE 

W.RTS 

REA  X IN 

SELECT 

I'JSt 

CPHS 

Driver 

EQLBR* 

FROZEN 

hcalc 

matrix 

ODES 

OUTi 

RKTOUT 

j 

i 

ROCKET 

1 

SAVE 

search 

SELECT 

K1NFU 

DRIVER 

frozen 

MAINJD 

ODES 

ODKINP 

OUTPUT 
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COMMON 

REFERRED  TO 

BLOCK 

BY  SUBROUTINE 

- 

OUT! 

PACK 

PRES 

REACT 

REAXIn 

RKTOUT 

ROCKET 

SELECT 

IKtM 

0ALIS 

MAlNj 

N2CEOM 

PICKUP- 

PROBLM 

8UMRV 

LKEQKN 

DRIVER 

ODES 

OOK 

PICKUP 

PROBLH 

RKTOUT 

ROCKET 

SURRY 

UKGEOH 

NZGEOM 

ODKINP 

SUHRV 

T02P 

LKMELT 

DERI  V 

BAIN1D 

LKO^K 

PICKUP 

SUMRY 

LKT8L 

PROBIM 

READS! 

LKTD2P 

AGP 

NEGEOH 

• 

PICKUP 

PROBL* 

REAOSI 

T02P 

LOOP 

BALI  8 

BLKOAI 

C0NTR2 

ERROR 

PIG! 
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COMMON 

BLOCK 


REFERRED  TO 
BY  SUBROUTINE 


COMMON 

BLOCK 

REFERRED  TO 
BY  SUBROUTINE 

ROCKET 

Save 

SEARCH 

MISCX 

CPHS 

EQLBRM 

FROZEN 

HCALC 

MATRIX 

OOES 

REACT 

RKTOUT 

ROCKET 

SAVE 

MOUHTS 

OUT  1 

ROCKET 

SEARCH 

MUST 

OERIV 

FLU 

NAMD 

AGP 

EISP 

!8P1D 

TD2P 

NAME* 

A0CALC 

CCALC 

^CALC 

JAMES 

PAR  TIL 

PCALC 

PROP 

TD2P 

TRACE 
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COMMON 

BLOCK 

REFERRED  TO 
BY  SUBROUTINE 

NAHfp 

ACOMP 

AGP 

AXI9PT 

CNTRL 

CON3T3 

EFN 

EISP 

I3P1D 

KPBPT 

ON  tO 

PARTlL 

PRINT 

PTINT 

3UMPJ 

SUMP2 

TBUDG 

T02R 

TRACE 

wall 

«DG! 

WLPT 

NAMEr 

ACO*P 

ONED 

PARTIL 

TD2P 

TRACE 

NAMEI 

ADJK 

CNTRL 

KP9PT 

PRINT 

T02P 

NAHEL 

AGP 

ONED 

PARTI L 

PROP 

TD2P 

TRACE 

wDGI 
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COMMON 

BLOCK 

REFERRED  TO 
BY  SUBROUTINE 

NAHEh 

ACOHP 

AXI3PT 

CHECK 

CNTRL 

COASTS 

KPBPT 

ONID 

PARTIt 

PRINT 

PTINT 

8UMP1 

8UMP2 

T02P 

TRACE 

WALt 

WDGI 

WLPT 

NAmLN 

ACOMP 

AxISPT 

EEN 

kpbpt 

PTINT 

SUMPt 

SUMP2 

T02P 

NAMEP 

ACOMP 

ADJK 

AXI3PT 

CHECK 

CNTPL 

COASTS 

EEN 

KPBPT 

PRINT 

PTINT 

SUMP  J 

8UHP2 

T02P 

Wl.PT 
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COMMON 

REFERRED  TO 

BLOCK 

BY  SUBROUTINE 

NAPL’Q 

ACOMP 

CNTRL 

COASTS 

ON  tO 

PARTIL 

PRINT 

Toap 

TRACE 

NAHrft 

ONtO 

partil 

PROP 

TD2P 

TRACE 

NAMES 

ONED 

PARTIL 

PROP 

Toap 

TRACE 

WOGI 

NAMET 

AXI3PT 

KPBPT 

PTINT 

SUMP! 

SUMP? 

TD2P 

NIPT 

NAMEW 

QDKINP* 

Tl'ZP 

WALL 

NAME* 

ONED 

PARTIL 

PRINT 

PROP 

T8LEDG 

TD2P 

TRACE 

nahev 

ONED 

PARTIL 

PROP 

TD2P 

TRACF 

NAMEZ 

PARTIL 

Toap 

TRACE 

COMMON 

BLOCK 


REFERRED  TO 
BY  SUBROUTINE 


NAME! 

NAME2 

NAME3 

NAME'* 

NAMIN 

NIN 


ABC ALC 

CCALC 

rcAtc 

JAMES 

ONED 

PARTU 

PCALC 

TD2P 

ABCALC 

CCALC 

ecalc 

JAMES 
PC  ALC 

ABCALC 

CCALC 

PCALC 

JAMES 

PARTIL 

PCALC 

PROP 

TRACE 

ABCALC 

CCALC 

MAIN1D  ** 
ODKINP* 

BLKDAT 

ERROR 

REAOIN 

STORE 


ODECUT 


OUT  1 

RKTOUT 

ROCKET 
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COMMON 

BLOCK 

REFERRED  TO 
BY  SUBROUTINE 

PERP 

DRIVER 

EQLBRM 

PROZEN 

ODES 

OJT 1 

RXTOUT 

ROCKET 

PEPPX 

iqlbrm 

PROZEN 

ODES 

OUT  1 

RXTOUT 

ROCKET 

PMPAN 

BAR3ET 

POINTS 

DRtvER 

EQLBRM 

PROZEN 

HCALC 

MATRIX 

ODES 

OUTJ 

RKTOUT 

ROCKET 

POINT* 

EQLBRM 

PROZEN 

HCALC 

MATRIX 

ODES 

RKTOUT 

ROCKET 

PTABLE 

PLU 

I AUX 

MAINtD 

ODK 

OCKINP 

PRES 

PT3AVE 

DRIVER 

DOES 

RDLARO 


DRIVER 

PROBtM 
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COMMON 

BLOCK 

RDJRF.X 

RVACO* 

SAVE 

3AV00E 

SHGAMA 

SPFCE3 

SPECfX 


REFERRED  TO 
BY  SUBROUTINE 


DRIVER 
ODES 
PACK 
REAX  JN 
SEARCH 

select 

STOICC 

INT 

MAINID 

ODK 

DDK  I HP 
OUTPUT 

CONTRJ 

FIG2 

PRXSM1 

DRIVER 

ODES 

DERI  V 
ODK 

ODK IMP 

CPHS 

EOLBRM 

P ROZEN 

HCAIC 

MATRIX 

ODES 

OUT  1 

RK  TOUT 

ROCKET 

SAVE 

SEARCH 


CPHS 

EOLBRM 

PROZEN 

HCALC 

MATRIX 

ODES 

OUT  l 

RK  TOUT 

ROCKET 

SA^E 

SEARCH 


COMMON 

REFERRED  TO 

BLOCK 

BY  SUBROUTINE 

3PINP0 

CPH3 

DRIVER 

EQLBRR 

ODES 

ODKINP 

OUT  1 

RACK 

READS! 

RKTOUT 

ROCKET 

search 

SELECT 

8UMRY 

TD2P 

TBLDTi 

DRIVER 
IAUX 
MAIN  JO 
OOKINP 

TBL0T2 

DRIVER 

IAUX 

HAINID 

TBLRC 

PICKUP 

SUMRV 

TBRSAV 

ODKINP 
Rt AX  IN 
8ELECT 

T020UT 

DRIVER 

PRINT 

TBIEDG 

T02P 

TD2RE 

PICKUP 

SUHRV 

TfcQfcC 

ODES 

THI3P 

PICKUP 

SUHRY 

tmwuat 

NZGEOH 

PICKUP 

8‘JMRY 

TD2P 

3-21 


COMMON 

REFERRED  TO 

BLOCK 

BY  SUBROUTINE 

T0T3C 

DRIVER 

IAUX 

mainid 

TRICON 

DRIVER 

8UMRY 

TRAJEC 

TSTA01 

DRIVER 

IAUX 

MAINJD 

TWOPHZ 

CONVRT 

FLU 

MAINID 

ODK 

ODKINR 

OUTPUT 

pack 

VERSON 

sumry 

WHERE 

CONTR1 

CUNTR2 

FIG1 

EIG2 

PRTSM1 

PRISM* 

WROOK 

ODK 

OUTPUT 

wtelow 

BARPRO 

DIRECT 

READSI 

ZDEBUG 

DEM  IV 
STF 

ZTR*N 

DRIVER 

IAUX 

MAINID 
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COMMON  BLOCK:  BALCOM 


FORTRAN 

NAME 

TYPE 

DIMENSION 

Description 

TIM 

L 

150 

Time,  tp  point  from  BAL  module 

PTC 

R 

150 

total  pressure  at  the  end  of  the  grain, 
Pci,  (PSIA) 

RSTAR 

R 

150 

r*  (ft) 

XLST 

R 

150 

L*  (In) 

ETACS 

R 

150 

*0*1 

XMDOT 

R 

150 

rrij,  mass  flow  rate  (slugs/sec) 

EMIEMP 

R 

150 

not  used 

PCI 

R 

150 

not  used 

XMT0T 

R 

150 

not  used 

FID 

R 

150 

not  used 

PMAX 

R 

- 

n.aximum  chambei  pressure  (PSIA) 

NBPT 

I 

- 

number  of  ballistics  time  points 

PAVE 

R 

- 

average  Pc  (PSIA) 

RMIN 

R 

- 

smallest  r*  in  the  RSTAR  array  (ft) 

XLAVE 

R 

- 

average  L*  (in) 

ECAVE 

R 

- 

average  r,c* 

XMAVE 

R 

- 

average  m (slugs/sec) 

Table  3-3  Common  Block:  BALC0M 
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COMMON  EEE 


F0RTRAN 

NAME 

TYPE 

DIMENSION 

DESCRIPH0N 

R 

15 

calculated  ^ 

R 

15 

empirical  rf 

ETAI 

R 

15 

Input  tj 

ETAU 

R 

15 

value  of  v to  be  used 

IPETA 

I 

15 

pointer  to  value  of  ^to  be  used 

Table  3-4  Common  Block  EEF. 
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C0MM0N  BL0CK: 

: LKEQKN 

F0RTRAN 

NAME 

DIMENSION 

DESCRIPTION 

IPUISP 

I 

- 

If  IPUISP  =1;  store  € and  ISP  In 
AEQR  and  XIEQR 

IREQ 

I 

- 

If  IREQ=1;  restricted  equilibrium  op- 
tion; if  0,  regular  equilibrium 

AEQR 

R 

50 

< 

^ corresponding  to  Isp  (XIEQR)  for 
the  restricted  equilibrliim  option 

XIEQR 

R 

50 

NEQR 

I 

- 

number  of  restricted  equilibrium  points 

SOLPT 

L 

- 

if  TRUE  solidification  point  proper- 
ties calculated 

FRZISP 

R 

- 

value  of  frozen  ISP 

Table  3-5  Common  Block  LKEQKN 
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COMMON  BLOCK:  LK0DK 

FORTRAN 

NAME 

TYPE 

DIMENSION 

Description 

AKIN 

R 

40 

^ corresponding  to  XIK(I) 

XIK 

R 

40 

kinetic  I calculated  by  ODK 

SD 

ETAKIN 

R 

- 

rtYihj  at  maximum  area  ratio  specified 

irnSe  $GEOM  namelist 

NKPT 

I 

- 

number  of  points  in  the  AKIN,  XTK  arrays 

Table  3-6  Common  Block:  LK0DK 


r 

I 

i 


[ 
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C0MM0N  BL0CK:  LKTBL 


F0RTRAN 

TYPE 

DESCRIPTION 

NAME 

IF0DE 

I 

If  IF0DE=1,  ODE  data  Is  available  to  the 
TBL  module 

IFTD2P 

I 

If  IFTD2P=1,  TD2P  data  is  available  to 

the  TD2P  module 

INUN 

I 

logical  unit  on  which  the  TD2P  data  Is 
available 

Table  3-  7 Common  Block  LKTBL 


1 


I 
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C0MM0N  HL0CK:  MAN 


J 

E 

l 


f 

I 


F0RTRAN 

DESCRIPH0N 

NAME 

0DF 

f if  value  = 0.0  not  considered 

BAL 

1 if  value  =1.0  module  by  that  name  is  to 

0DK 

J be  executed 

TD2P 

^ if  value  = 2.0  module  data  was  previously 

j generated 

TBL 

V if  value  = 3.0  loss  to  be  read  in 

Table  3-8  Common  Block  MAN 


| 


l 

I 
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C0MM0N  HL0CK:  MARKET 


F0RTRAN 

NAME 

TYPE 

DESCRIPTION 

LODE 

L 

If  .TRUE,  module 

ODE 

is  to  be  executed 

LODK 

L 

if  .TRUE,  module 

ODK 

is  to  be  executed 

LTD2P 

L 

if  .TRUE,  module 

TD2P 

is  to  be  executed 

LTBL 

L 

if  .TRUE,  module 

TBL 

is  to  be  executed 

LBAL 

L 

if  .TRUE,  module 

3AL 

is  to  be  executed 

Table  3-9  Common  Block  MARKET 


I' 

l 
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C0MM0N  BL0C: 

TBLRE 

FORTRAN 

NAME 

TYPE 

DIMENSION 

DESCRIPTION 

NTBL 

l 

- 

number  of  TBL  points 

ABL 

R 

150 

A/A*  from  TBL 

DELISP 

R 

150 

A spTBL 

DISP 

R 

^1  at  the  maximum  area 

spTBL 

ratio  specified  In  the  $GEOM 
namelist 

Table  3-10  Common  Block:  TBLRE 


i 
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COMMON  BLOCK:  TD2RE 


FORTRAN 

NAME 

TYPE 

DIMENSION 

Description 

A2D 

R 

150 

q corresponding  to  values  of  XI2D  and  E2D 

XI2D 

R 

150 

I ((.)  from  TD2P 
* 

E2D 

R 

150 

r?TD2P  W 

NTDPT 

I 

number  of  points  in  the  A2D,  XI2D  and  E2D 
array 

ETD2P 

R 

- 

time  average  value  of  t)td2P 

XITD2P 

R 

- 

time  averaged  value  of  I 

spTD2P 

ETD2PA 

R 

- 

7lTD2P  corrected  for  ^roat  erosion 

CD 

R 

- 

nozzle  discharge  coefficient  as  calculated 
by  TD2P 

WD0T2D 

R 

nozzle  mass  flow  rate  as  calculated  by 
TD2P 

Table  3-11  Common  Block:  TD2RE 
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COMMON  THISP 


F0RTRAN 

TYPE 

DESCRIPTION 

NAME 

ISP 

R 

Theoretical  I at  maximum  area  ratio 
specified  In  me  $GEOM  namelist 

CSTAR 

R 

Theoretical  c*  (ft/sec)  as  calculated  by 
the  0DE  module 

Table  3-12  Common  Block  THISP 
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COMMON  HLOCK:  THROAT 


FORTRAN 

NAME 

TYPE 

DESCRIPTION 

DRSDT 

R 

dr*/dt  (ft/sec) 

TBIJRN 

R 

bum  time  (sec) 

NAR 

I 

number  of  points,  at  which  nozzle  erosion  data  is  known 

FSUBM 

! R 

length  of  submergence/length  of  internal  motor 

AEAS 

R 

nozzle  entrance  area/nozzle  throat  area 

AVELS 

R 

motor  average  L*  (in) 

KN0S 

I 

nozzle  type  flag;  =1  steel  nozzle,  =0  all  others 

FCONP 

R 

, mole  fraction  of  condensed  phase  (moles/100  grs  of  product) 

COMMON  BLOCK:  TRJCOM 

FORTRAN 

NAME 

TYPE 

DIMENSION 

Description 

ITJF 

I 

- 

trajectory  flag,  if  0 no  trajectory  input 

PAMBI 

R 

51 

ambient  pressure  (PSIA) 

TTRJ 

R 

51 

times  corresponding  to  PAMBI 

NTJPT 

I 

- 

number  of  trajectory  points  input 

Table  3-14  Common  Block  TRJCOM 


r 

f. 


I 
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4.0  Program  Files 


The  SPP  program  uses  nine  files  to  read,  write,  punch  and  temporarily  save 
data.  Six  of  these  files  are  always  referred  to  by  variable  names  and  hence,  may 
be  changed  to  suit  the  restrictions  of  the  computer  system  being  used.  These 
variable  names  are  set  in  Subroutine  DRIVER  for  units  25  through  29  and  by  input 
data  for  unit  3 in  subroutine  NZGE0M. 

Logical  unit  8 is  equated  to  the  system  punch  file  for  CDC  6000  series 
computers  and  it  is  rewound.  Therefore,  care  must  be  taken  in  assigning  this 
file  on  machines  which  do  not  allow  rewinds  on  the  system  punch  file  (e.g.,  Univac 
1108).  This  linkage  data  written  out  on  this  unit  is  described  in  Volume  III,  Section 
2.11  and  will  not  be  repeated  here. 

Table  4-1  lists  the  files  used  by  the  SPP  code  and  the  FORTRAN  variable 
names  used  to  reference  these  files.  The  device  types  listed  correspond  with  the 
current  usage  on  CDC  6000  series  computers.  The  mnemoics  used  in  Table  4-1 
are  listed  below. 

T - Tape  files 

MS  - Mass  Storage  (perm)  files 

TMS  - Temporary  Mass  Storage  (scratch)  files 
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FORTRAN 

VARIABLE 

NAME 

Loc^cal 

TJnit 

Device  Type 

Description 

INP 

3 

T,  MS,  TMS 

alternate  Input  file  for  SPP  linkage 
data 

- 

5 

card  reader 

data  input  sLream 

- 

6 

line  printer 

printed  output  stream 

- 

8 

punch 

module  linkage  data  output  file 

JANAF 

25 

T,  MS,  TMS 

thermodynamic  data  file 

KREAX 

2f> 

TMS 

temporary  file  used  to  store  reaction 
rate  data 

KSTF 

27 

TMS 

temporary  file  used  to  store  a short 
version  of  the  thermodynamic  data 

IRREAD 

28 

TMS 

temporary  file  used  for  multiple  cases 
in  tne  0DK  module. 

ITSTAB 

29 

TMS 

the  data  written  on  this  file  is  not 
currently  used 

Table  4-1  Files  U3ed  By  The  SPP  Code 


i 
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5.0  PROGRAM  SUBROUTINE  DESCRIPTIONS 


This  section  of  Volume  II  contains  a subroutine  by  subroutine  description 
of  the  SPP  computer  code.  Dictionaries  of  some  of  the  more  important  variables 
which  appear  in  labeled  common  areas  of  the  program  appear  in  Section  3. 

Section  5 is  organized  in  accordance  with  the  overlay  structure  of  the  SPP 
code  and  the  user  is  advised  to  refer  to  Section  3 of  this  document  for  the  location 
of  individual  subroutines  and  common  blocks  and  their  cross  references. 
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5.1  Link  00  Subroutines 


5.1,1  Program  Driver 

This  is  the  main  routine  for  the  SPP  code  and  as  such  provides  the  overlay 
control , defines  most  of  the  upper  level  labeled  common  blocks,  and  initializes 
most  of  the  control  variables  used  in  the  code.  On  CDC  6000  computers  this 
routine  also  defines  all  of  the  files  and  buffer  sizes  which  are  used  in  this  code. 

The  DRIVER  routine  is  essentially  driven  by  directive  cards  which  are  read 
in  on  logical  unit  5.  The  directive  cards  which  are  recognized  by  this  subprogram 
are  listed  below  in  Table  5.1-1  along  with  the  action  taken  by  the  DRIVER  routine. 


Table  5.1-1  Data  Input  Directives 


Directive  card 

Action  Taken 

THERM  0 

Call  in  overlay  1,0  which  processes  ther- 
modynamic data. 

TITLE 

Stores  the  case  title  for  subsequent  print 
out. 

PR0BLEM 

Cells  in  the  PR0BLM  subroutine  which  de- 
termines which  loss  modules  are  to  be  ex- 
ecuted. This  directive  and  subsequent 
data  initiate  the  execution  of  the  various 
loss  modules. 

LOW  T CPHS 

Calls  in  subroutine  LTCPHS  which  processes 
low  temperature  thermodynamic  data. 

GE0METRY 

Calls  in  subroutine  GE0M  yvhich  processes 
the  nozzle  geometry  data. 

TRAJECT0RY 

Call  in  subroutine  TRAJ  which  processes 
the  trajectory  and/or  ambient  pressure 
history  of  the  motor. 

In  addition  to  processing  the  directive  cards  and  calling  in  the  selected 
loss  modules,  the  main  program  also  hi  idles  the  communication  of  the  linkage 
data  between  modules  and  calls  the  summary  data  output  routine  SUMRY. 
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5.  1 . 2 Subroutine  BL0CK  DATA 


BL0CK  DATA  contains  atomic  data  stored  in  AT0M(i,j)  and  many  of  the  vari- 
ables used  with  the  variable  format,  FMT.  The  AT0M  variables  are  defined  in 
appendix  B,  Reference  1.  The  format  variables  are  stored  in  the  common  labeled 
0UPT  and  are  described  here. 

A variable  format  was  used  so  that  one  format,  FMT,  could  be  used  in  the 
final  output  with  changes  in  the  number  of  decimal  places  according  to  the  sizes 
of  the  numbers.  The  format  is  used  to  print  a label  and  from  1 to  13  associated 
numbers.  The  labels  contain  14  alphameric  characters  stored  in  four  words  and 
printed  with  3A4,  A2.  The  numbers  are  all  printed  in  a field  of  9.  FMT  is  initially 
set  in  HL0CK  DATA  as  follows: 


FMT 

(i) 

(2) 

(3) 

(4) 

(5) 

(6) 

(7) 

(8)  (9) 

do) 

(11) 

(1H 

,3A4 

,A2, 

F9. 

0, 

F9. 

0, 

F9.  0.. 

F9. 

0, 

FMT 

(12) 

(13) 

(14) 

(15) 

(16) 

(17) 

(18) 

(19) 

(20) 

(21) 

F9. 

0, 

F9. 

0, 

F9. 

0, 

F9. 

0, 

F9. 

0, 

FMT 

(22) 

(23) 

(24) 

(25) 

(26) 

(27) 

(28) 

(29) 

(30) 

F9. 

0, 

F9. 

0, 

F9. 

0, 

F9. 

0 

) . 

where  the  spaces  are  stored  as  blanks. 


Some  variables  set  in  BL0CK  DATA  to  modify  FMT  are  as  follows: 
Variable:  FO  FI  F2  F3  F4  F5  FB  FMT13  FMT9X  FMT19 

Storage:  0,  1,  3,  4,  5,  13,  9X,  19, 

The  following  is  a list  of  variables  used  as  labels  and  printed  with  3A4, 

A2  in  FMT: 


Variable 

Stored  label 

FP 

P,  ATM 

FT 

T,  DEG  K 

FH 

H,  CAL  G 

FS 

S,  CAL/(G)  (K) 

FM 

M,  MOL  WT 

FV 

(DLV/DLP)  T 

FD 

(DLV/DLT)  P 

FC 

CP,  CAL/ (G)  (K) 

FG 

GAMMA  (S) 

FL 

SONVEL,  M/SEC 

FR1 

PC/P 

FC1 

CF 

FN 

MACH  NUMBER 

FR 

CSTAR,  FT/SEC 

FI 

ISP,  LB-SEC/LB 

FA 

IVAC,  lb-secab 

FA1 , FA2 

AE/AT 
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5.1.3  Suoroutine  ATM0S 


Soubroutine  ATM0S  is  used  by  subroutine  TRAJEC  (see  Section  5.1.16)  to 
calculate  the  ambient  pressure  history  during  a motor  firing  if  the  ambient  condi- 
tion time  history  is  given  as  a function  of  altitude. 

This  subprogram  calculates,  given  the  geometric  altitude,  the  ambient  temp- 
erature PR),  pressure  (psia),  sound  speed  (ft/sec),  density  (slugs/ft3),  and  pres- 
sure gradient  (lb/ft3).  Subroutine  ATM03  assumes  that  the  atmosphere  consists 
of  a gas  of  constant  molecular  weight  on  a spherical  earth  whose  atmosphere  is 
in  a hydrostatic  balance.  The  temperature  profiles  used  in  the  hydrostatic  balance 
are  those  given  in  the  1962  ARDC  standard  atmosphere  tables. 

These  assumptions  yield  data  which  agree  with  that  standard  atmosphere 
table  to  within  1%  up  to  altitudes  of  700  kilometers  (12,296,585  ft.). 


5.1.4  Subroutine  AVEBAL 

Subroutine  AVE8AL  is  called  by  subroutine  PICKUP  and  calculates  from  the 

linkage  data  supplied  by  the  internal  ballistics  calculation  the  time  averaged  values 

of  chamber  pressure,  L*,  throat  radius,  rj  *,  and  mass  flow  rate.  Also  printed  out 

c 

are  the  maximum  and  minimum  values  of  the  above  quantities  along  with  their  time 
histories.  The  throat  erosion  rate  in  inches/second  is  also  calculated  and  printed 
out. 
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5.1.5  Subroutine  ETACFX 


This  subroutine  calculates  the  nozzle  performance  losses  using  empirical 
correlations.  The  losses  calculated  in  this  routine  are  in  terms  of  per-cent  nozzle 
thrust  coefficient  loss.  Subroutine  SUMRY  converts  these  losses  into  fractional 
Igp  efficiencies  using  the  equation 

100^ 

*ISV1  = Too  *»c*  5,1_1 

where  ^is  the  nozzle  Cp  loss  as  calculated  below 

Care  should  be  taken  when  using  the  above  equation  that  c*  efficiency, 

n *,  is  not  used  more  than  once  to  convert  Cr  efficiencies  to  I efficiencies. 

"C1*  r Sp 

The  losses  considered  and  equations  used  by  this  routine  are  given  below. 

A more  complete  description  of  these  losses  is  given  in  Voi.  I,  Section  4. 


Nozzle  divergence  loss  is  calculated  as: 


5.1-2 


where  <y=half  angle  of  conical  nozzle,  included  angle  of  contoured  nozzle 
0Ex  =exit  angle  of  contoured  nozzle 

The  loss  due  to  finite  rate  chemical  kinetics  is  calculated  as: 


^KIN  = 33,3 


theoretical  frozen  Isp  ,200  \ 
theoretical  shifting  Isp  * P 


where  P=nczzle  chamber  pressure  in  psia  or  200,  whichever  is  greater. 


The  boundary  layer  loss  is  expressed  in  the  form  of  a time-dependent  heat 
transfer  expression,  with  a correction  term  for  nozzle  expansion  ratio: 


,0.8 


^BL  " C1 


D 


0.2 


1 + 2 exp  (-C2P0,8t/D°'2) 


1 + 0.016U-  9) 


5.1-4 


where  P = pressure,  psi 

= throat  diameter,  in 
t - time,  sec 
e - expansion  ratio 
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The  coefficients  and  C2are  determined  as  Is  shown  below.  The  nozzle 
type  is  input  by  the  user  as  KN0Z  in  the  $GE0M  namelist. 

Ordinary  Nozzles:  = 0.00365 

C2  = 0.000937 

Steel  Nozzles:  = 0.00506 


The  nozzle  two  phase  flow  loss  term  is  calculated  as: 

/4Dp°5 

nTP  p0.15  f0.08D  Cg 

where  $ = mole  fraction  of  condensed  phase,  moles/100  gm 

Dp=  particle  size,  u 
P = pressure,  psi 
<r  •=  expansion  ratio 
= throat  diameter,  in. 


The  coefficients  are  determined  from  the  following: 
If  £ ^ 0.09,  C4  = 0.5 


with  Dp  =•  0.454  P1//3  £1//3  [l  - exp  (-0.004L*)  J [ 1 + 
and  L*  = motor  characteristic  length,  in. 

The  motor  submergence  loss ‘is  calculated  by: 

pi,  0.8  qO  • 4 

"SUB  = 0*0684  (A^  XT 

Ut 

where 

P = pressure,  psi 

£ = mole  fraction  of  condensed  phase,  moles/100  gm 

A*  = nozzle  entrance  area/nozzle  throat  area 
S = length  of  submergence/length  of  internal  motor 
Dt  = throat  diameter,  in. 


I 

I 


i 

f 

t 


0.045DtJ  5.1-6 


5.1-7 
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5,1,6  Subroutine  FIND 


This  subroutine  locates  the  indes,  I,  in  a table  such  that  X(I)  <X  <X(I+1). 
On  option,  when  the  variable  IEXTP  in  labeled  common  IXTRAP^O,  subroutine  FIND 
allows  for  extrapolation  in  tables.  See  Appendix  A for  conversion  aides. 


5,1.7  Subroutine  LINK1 A 

Subroutine  LINKlA's  only  function  is  to  call  subroutine  PR0BLM.  It  is 
provided  to  allow  addition  overlay  capability  on  computers  with  explicid  overlay 
calling  features. 


5.1.8  Subroutine  L7CPHS 

This  subroutine  processes  the  low  temperature  C , H,  S Thermodynamic 
Data  extension  input  as  described  in  detail  in  Volume  III,  Section  2.2.1.  The 
low  temperature  values  of  Cp,  H,  and  S are  at  present  only  used  in  the  0DK 
module. 
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5,1,9  Subroutine  NZGE0M 

Subroutine  NZGE0M  Is  called  by  program  DRIVER  when  a GE0METRY  directive 
is  encountered  In  the  card  input  stream.  This  routine  uses  the  $GE0M  Namelist 
input  to  define  the  nozzle  geometry.  Data  which  are  read  under  this  Namelist 
are  stored  in  labeled  common  blocks  LKGE0M,  LKBM,  ARPRNT,  BALC0M,  LKTD2P, 
and  THR0AT. 

The  variables  in  the  above  common  blocks  are  described  in  Section  3,2, 
However,  for  completeness,  the  variable  names  whose  values  are  determined  in 
subroutine  NZGE0M  are  listed  below  in  Table  5,1-1, 


C0MM0N  HL0CK 

VARIABLE  NAME 

ARPRNT 

ASUB 

AS  UP 

NASUB 

NAS  UP 

BALC0M  ~~ 

RSTAR 

TIM 

LKBM 

NN0Z  " 

R 

LKGE0M 

GMVAR 

IGM 

IWF 

NWS 

RS 

ZS 

LKTD2P 

RSTARM 

XLS 

THR0AT 

DRSDT 

FSUBM 

NAR 

TBURN 

Table  5,1-1  Common  Variables  Set  in  NZGE0M 

The  variables  shown  in  the  above  table  are  communicated  downward  only 
in  the  SPP  code.  That  is,  values  set  in  this  routine  may  be  changed  in  indivi- 
dual calculational  mcdules,  however,  the  values  which  are  changed  are  not 
transmitted  to  other  calculation  modules.  Instead  the  initial  values  set  via  the 
$GE0M  Namelist  are  used. 
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5.1,10  Subroutine  PICKUP 


Subroutine  PICKUP  processes  all  of  the  linkage  data  between  the  five  basic 
computational  modules  which  are  transmitted  on  logical  unit  8 except  for  the 
boundary  layer  edge  conditions  generated  by  the  TD2P  module  for  the  TBL  module. 

This  routine  also  performs  unit  conversion  on  the  data  read  in,  checks  for 
overflow  of  tables,  and  prints  out  a summary  of  the  results  for  each  of  the  modules 
which  have  been  executed. 


5.1.11  Subroutine  FR0BLM 

Subroutine  PR0BLM  is  executed  by  program  DRIVER  when  the  PROBLEM 
directive  is  encountered  in  the  card  input  stream. 

This  routine  sets  flags  and  determines  via  user  input  in  the  $PR0B  Namelist 
which  of  the  five  basic  computational  loss  modules  are  to  be  executed.  If  the 
option  to  read  previously  generated  data  is  selected,  subroutine  PICKUP  is  called 
by  PR0BLM  to  read  these  data. 


5.1.12  Subroutine  SETUP 

Subroutine  SETUP  is  used  to  position  the  linkage  data  file,  logical  unit  8, 
for  processing  by  subroutine  PICKUP.  This  routine  reads  through  the  linkage  data 
file  until  the  appropriate  calculation  module  name  is  found,  then  it  backspaces 
one  logical  record  so  that  the  loss  module  name  is  the  first  card  image  read  by 
subroutine  PICKUP. 
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5.1.13  Subroutine  SPLN 


Performs  either  cubic  or  linear  Interpolation  between  two  given  points. 


Cubic  interpolation  for  a function  and  Its  first  two  derivatives  Is  performed 
as  described  below:  Given  function  values  yn  and  yfi+1  and  first  derivative  values 
y'n  and  y^+1  at  xR  and  xn+1,  this  subroutine  evaluates  y(x),  y'(x),  and  y"(x)  for 
xn  <:x<xn+1  using: 

y = A(x  - x J3  + B(x  - x )s  + C(x  - x ) + D 
n n n 


V-  (y;H  - v')  A 

where  - 

A = F • [^n+l+^h-*J 

B =”F  * [ 'yn+l  + 2yn)h  ' 3k] 

c y; 

D=  yn 
h = Vl  • xn 


Linear  interpolation  for  a function  and  its  first  two  derivatives  is  performed 
as  described  below: 


V' 


Vl  I yn 

xn+l  " Xn 


y"  = 0.0 
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5.1,14  Subroutine  SUMRY 


This  subroutine  is  called  by  program  DRIVER  (Section  5.1.1)  and  prints  out 
a summary  report  on  the  results  of  each  of  the  five  basic  computational  modules 
(if  executed)  at  the  maximum  area  ratio  requested  in  the  $GE0M  namelist.  If 
the  GE0METRY  directive  was  not  input  to  the  SPP  code  the  SUMRY  routine  assumes 
that  one  or  more  of  rhe  calculation  modules  were  executed  independently  and 
that  no  summary  of  results  is  desired. 

There  are  three  basic  sources  for  each  of  the  losses  considered  by  the  SPP 
code.  These  sources  are  due  to  one  of  the  following  items  shown  in  the  table  below. 


1.  User  Input 

2.  Execution  of  a loss  module 

3.  Calculation  of  an  empirical  correlation 


Table  5.1-3  Sources  of  Losses  Considered 

User  input  of  loss  efficiencies  along  with  selection  of  which  of  the  loss 
modules  are  to  be  executed  is  determined  by  input  to  subroutine  PR0BLM  (Section 
5.1.11).  Empirical  losses  are  calculated  by  subroutine  ETACFX  (Section  5.1.5) 
which  is  called  by  SUMRY.  While  all  of  the  above  losses  are  printed  out  in  this 
routine,  the  selection  of  which  source  of  the  loss  is  to  be  used  in  determining  the 
total  delivered  perfoimance  of  the  solid  propellant  rocket  motor  is  given  in  the 
order  appearing  in  Table  5.1-3.  That  is,  if  the  user  of  the  SPP  code  inputs  a loss, 
that  is  the  loss  efficiency  that  is  used  in  calculating  the  final  delivered  perfor- 
mance. Hence,  while  user  input  overrides  values  calculated  by  the  analytical 
loss  modules,  the  values  calculated  by  the  loss  modules  takes  precedence  over  the 
values  determined  by  the  empirical  correlations. 

Using  the  selection  criteria  described  above,  subroutine  SUMRY  then  prints 

out  the  product  of  efficiencies  along  with  the  decrement  of  performance  due  to  a 

turbulent  boundary  layer.  Hence,  if  no  input  efficiencies  were  input,  the  value 

of  the  delivered  performance  due  to  input  losses  would  then  be  just  the  theoretical 

I calculated  by  the  0DE  module, 
sp 


In  order  to  determine  the  motor  delivered  performance  to  ambient  conditions, 
trapezodial  rule  integration  is  used  to  calculate  the  average  ambient  pressure  and 
average  expansion  ratio  if  this  information  is  available.  That  is,  the  average 
performance  delivered  to  ambient  conditions  is  calculated  by 

_ n 

I =1  f),  . - &I  - A I 5.1-8 

SpD  spth  1=1  ‘l'k  SpTBL  spA 

where  I = I calculated  by  0DE 

spth  spVAC 

v . = fractional  loss  efficiency  determined  via  the  selection  rule, 

k,  for  the  Ith  loss  considered 

A I = loss  decrement  due  to  the  formation  of  a turbulent  boundary 
spTLL  layer 

A I = loss  decrement  due  to  the  motor  delivering  thrust  to  an  am- 
spA  blent  pressure.  This  value  is  calculated  as  (P  . /?) 
x 7 xcvlg^Cj^ 

If  both  the  ballistics  and  the  two  dimensional -two  phase  flow  modules  link- 
age data  are  available,  then  subroutine  SUMRY  calculates  the  delivered  perfor- 
mance throughout  the  motor  firing  for  each  time  point  specified  in  the  trajectory 
and  ballistics  input  data.  The  following  formulas  are  used  to  compute  the  values 
printed  out  by  this  routine. 

nCE(t)  = max  (1.0,  CD  * rrt.*(t)) 

spTD2P^ 

IspVAC(t)  = Ispth  * f,®2pW  ' 4spTD2p(t=0))  * ’7CE<t)  ‘ ”SUB  * nKIN  " ^TBL 


FVAC(t) 

= i * mft) 

spVAC 

Famb(t) 

= Pamb^  - Aexlt 

Ispamb(y 

= Fa,b(t)/m(t) 

fdel 

= FVAC(t)  " Famb(t) 

IspDEL 

IspVAC  Ispamb 

I 

= J FDEL  ^ dt 
0 

m 

= J*  mdt 

0 
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where 


m 


nc* 


amb 


^SUB 

%IN 


spbl 


^TD2P 

SpTD2P 

A .* 
exit 

*amb 

I 


m 

I 

spDEL 

F 

DEL 

I 

spamb 


mass  flow  rate  (from  Ballistics  calculation) 
c*  efficiency  (from  Ballistics  calculation) 
ambient  pressure  (from  input  tables) 

submergence  loss  efficiency  (selected  as  described  above) 
kinetic  loss  efficiency  (selected  as  described  above) 
boundary  layer  loss  decrement  (selected  as  described  above) 

2D-2phase  flow  loss  efficiency  (from  TD2P  as  a function  of  *(t)) 

2D-2phase  I (from  TD2P  as  a function  of  ^ (t)) 

sp 

no?zle  exit  area  (from  input) 

thrust  decrement  due  to  ambient  pressure 

total  impulse  to  this  time 
total  mass  expended  to  this  time 
delivered  specific  impulse 

delivered  thrust 

I decrement  due  to  ambient  pressure 
sp 
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5.1,15  Subroutine  TIMERX 


This  subroutine  Is  provided  to  localize  printing  of  computer  execution  times 
during  the  calculations.  This  subroutine  should  be  modified  at  each  local  instal- 
lation to  provide  proper  interface  with  local  system  timing  routines.  This  subroutine 
calls  a CDC  system  subroutine  named  SEC0ND  to  obtain  the  current  run  execution 
time  in  seconds. 


5.1.16  Subroutine  TRAJEC 

Subroutine  TRA/EC  is  used  to  calculate  ambient  pressure  during  a motor  firing. 
A table  of  ambient  pressure  or  altitude  is  input  to  this  routine  via  the  Namelist 
name  $TRAJ  as  a function  of  time.  Delivered  performance  to  the  specified  amiJent 
conditions  is  then  calculated  in  subroutine  SUMRY  and  printed.  If  an  altitude 
versus  time  table  is  input  to  subroutine  TRAJEC,  then  subroutine  ATM0S  is  used  to 
calculate  the  ambient  pressure  at  which  the  performance  of  the  motor  is  delivered. 
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5.2  Link  10  Subroutines 


5.2.1  Subroutine  LINK10 

This  subroutine  consists  only  of  a call  to  subroutine  TTAPE  and  is  used 
only  to  facilitate  conversion  of  the  program  overlay  structure  to  other  computers. 


5.2.2  Subroutine  TTAPE 

When  a THERM#  directive  card  is  read  by  the  main  program  this  subroutine 
is  called  to  generate  a master  Thermodynamic  Data  tape  (Logical  Unit  25).  The 
input  Thermodynamic  Data  is  in  curve  fit  form  and  is  identical  to  that  required  for 
the  #DE  computer  program  described  in  NASA  SP-273,  Reference  1.  The  format 
for  this  curve  fit  date  is  described  in  the  Volume  III,  Section  2.2.  This  routine 
does  not  update  the  master  Thermodynamic  Data  tape,  but  instead  creates  a new 
one.  Hence,  all  of  the  thermodynamic  must  be  read  in,  old  and  new,  if  the  user 
wishes  to  add  new  species  to  this  file. 
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5.3  Link  20  Subroutines 

Overlay  Link  20  contains  the  0DE  module  which  calculates  the  theoretical 

or  reference  I used  in  the  methodology  incorporated  in  the  SPP  code, 
sp 

The  main  body  of  subprograms  which  are  contained  in  this  overlay  are  from 
the  NASA-Lewis  Equilibrium  Program  (Ref.l)  as  modified  for  rocket  performance 
calculations  (Ref.  3).  The  brief  analytical  descriptions  of  the  equilibrium  calcu- 
lations contained  herein  were  condensed  from  Reference  1.  For  a complete  des- 
cription of  these  calculations  see  Reference  1 (NASA  SP-273). 


5.3.1  Subroutine  LINK20 

This  subroutine  consists  only  of  a call  to  subroutine  0DES  and  is  used  only 
to  facilitate  conversion  of  the  program  overlay  structure  to  other  computers. 
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man 


5.3.2  Subroutine  CPHS 


Cp° 


This  subroutine  evaluates  the  thermodynamic  functions 


HX 


S°T 


RT  ' R 

from  curve  fit  coefficients.  Two  sets  of  coefficients  are  used  for  two  adjacent 
temperature  Intervals.  The  functions  evaluated  are  presented  below: 

Cp\ 


~ 1T 

h°t 

~kT 

S T 
R 

G°T 

TT 


a,  + a0T  + a.T2  + a ,T3  + a.T4 
1 2 3 4 5 


ai  + 


a2T 


a3T" 


a1  InT  + a2T  + 


a3T" 


a.r 


a.T3 

4 


a5T4 


T 


a.T4 

+ "i4_  + a7 


H°T 

RT 


^T 

R 


5.3.3  Subroutine  DETON 

This  is  a dummy  subroutine  used  to  replace  subroutine  DET0N  of  Reference 

1. 


5.3 .4  Subroutine  EFMT 

Subroutine  EFMT  (E-format)  writes  statements  in  a special  exponent  form. 
This  form  is  similar  to  the  standard  FORTRAN  E-format,  but  the  letter  E and  some 
of  the  spaces  have  been  removed  for  compactness.  This  routine  is  not  used  in 
the  SPP  code  at  the  present  time. 
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5,3,5  Subroutine  EQLBRM 


EQLBRM  is  the  control  routine  for  the  equilibrium  module  which  calculates 
equilibrium  compositions  and  thermodynamic  properties  for  a particular  point.  A 
free-energy  minimization  technique  Is  used.  This  routine  permits  calculations 
such  as  (1)  chemical  equilibrium  for  assigned  thermodynamic  states  (T, P) , (H,P), 
(S,P),  (T,V),  (U,V),  or  (S,V),  (2)  theoretical  rocket  performance  for  both  equili- 
birum  and  frozen  compositions  during  expansion,  (3)  incident  and  reflected  shock 
properties,  and  (4)  Chapman- Jouguet  detonation  properties.  The  program  considers 
condensed  species  as  well  as  gaseous  species.  A detailed  description  of  the 
equations  and  computer  program  for  computations  involving  chemical  equilibria 
in  complex  systems  is  given  in  Reference  1.  Figures  4(a)  through  4(c)  of  Reference 
1 give  a complete  flow  diagram  for  this  subroutine. 


5,3,6  Subroutine  FR0ZEN 

Subroutine  FR0ZEN  is  called  from  subroutine  R0CKET  to  calculate  the 
temperature  and  thermodynamic  properties  for  the  following  assigned  conditions: 

1.  Composition  frozen  at  combustion  conditions. 

2.  At  assigned  exit  pressure. 

3.  At  assigned  entropy  equal  to  the  entropy  at  combustion  conditions. 

The  iteration  procedure  used  for  obtaining  the  exit  temperature  is  discussed  in  the 
section  Procedure  for  Obtaining  Frozen  Rocket  Performance  (p.  40,  Reference  1). 

This  subroutine  has  been  modified  to  calculate  the  strictly  frozen  performance 
at  well  below  the  solidification  point  of  any  condensed  phases  which  are  present. 


5-18 


5.3.7  Subroutine  GAUSS 


Subroutine  GAUSS  Is  used  to  solve  the  set  of  simultaneous  linear  Iteration 
equations  constructed  by  subroutine  MATRIX.  The  solution  Is  effected  by  perform- 
ing a Gauss  reduction  using  a modified  pivot  technique.  In  this  modified  pivot 
technique  only  rows  are  interchanged.  The  row  to  be  used  for  the  elimination  of 
a variable  is  selected  on  the  basis  that  the  largest  of  its  elements,  after  divi- 
sion by  the  leading  element,  must  be  smaller  than  the  largest  element  of  the  other 
rows  after  division  by  their  leading  elements. 

The  solution  vector  is  stored  In  X(k),  In  the  event  of  a singularity,  IMAT 
(which  is  equal  to  the  number  of  rows)  is  set  equal  to  IMAT  - 1.  IMAT  is  tested 
later  in  subroutine  EOLBRM. 


5 . 3 . 8 Subroutine  HCALC 


The  purpose  of  HCALC  is  to  calculate  thermodynamic  properties  for  reactants 
under  certain  circumstances.  HCALC  is  called  from  entry  NEW0F  of  SAVE. 


HCALC  is  called  from  NEW0F  when  CALCH  is  set  true . CALCH  is  set  true 
in  the  main  program  when  zeros  have  been  punched  in  card  columns  37  and  38  on 
one  or  more  REACTANTS  cards.  The  zeros  are  a code  indicating  that  the  enthalpy 
(or  internal  energy  for  UV  problems)  for  the  reactant  should  be  calculated  from  the 
THERM0  data  at  the  temperature  punched  on  the  card.  This  temperature  has  been 
stored  in  the  RTEMP  array.  CPHS  is  called  to  calculate  the  enthalpy.  The  value 
is  stored  in  the  ENTK  array  and  printed  in  the  final  tables. 

The  properties  calculated  in  subroutine  HCALC,  their  FORTRAN  symbols, 
and  the  conditior  for  which  they  are  used  are  as  follows: 


Property 

F0RTRAN  Symbol 

When  used 

H(K)T 

HPP(K) 

if  cc  37  and  38  are 

VR 

HSUBO 

input  as  00  on  reactants  card 

5.3.9  Subroutine  MATCH 

This  subroutine  is  called  by  R0CKET  to  supply  subroutine  MUK  with  a vector 
of  internal  sequence  m.mbers  which  point  to  the  appropriate  Lennard-Jones  parameters 
used  by  MUK  to  calculate  transport  properties.  The  input  to  MATCH  is  the  species 
name  (from  the  SUB  array)  and  the  output  corresponds  to  the  index  numbers  in  the 
table  in  the  MUK  write  up. 
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5.3.10  Subroutine  MATRIX 


This  subroutine-  sets  up  the  matrices  corresponding  to  tables  I through  IV 
of  Reference  1.  The  assigned  thermodynamic  state  being  set  up  (tables  I and  II) 
is  specified  by  the  following  codes: 


Assigned 

thermodynamic 

state 

Codes 

TP 

TP  = .TRUE.  VOL  = .FALSE.  CONVG  = .FALSE. 

HP 

HP=  .TRUE.  VOL  = .FALSE.  CONVG  = .FALSE. 

SP 

SP  = .TRUE.  VOL  = .FALSE.  CONVG  = .FALSE. 

TV 

TP  = .TRUE.  VOL  = .TRUE.  CONVG  = .FALSE. 

UV 

HP=  .TRUE.  VOL  = .TRUE.  CONVG  = .FALSE. 

SV 

SP  = .TRUE.  VOL  = .TRUE.  CONVG  = .FALSE. 

After  convergence  of  any  of  the  previous  six  problems,  setup  of  the  deriva- 
tive matrices  (tables  III  and  IV)  is  specified  by  the  following  codes: 


Derivative 

Codes 

DLVTP 

DLVPT 

1 1 

Note:  Only  the  Rocket  option  is  available  in  the  SPP  code. 
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5.3.11  Subroutine  MUK  (CP,  XW,  T,  C,  IT,  N,  XM,  XMT,  TK,  PR) 

This  routine  calculates  the  viscosity,  thermal  conductivity  and  Prandtl  number 
for  a gas  composed  of  a mixture  of  species. 

CALLING  SEQUENCE; 


CP 

is  an  array  of  N specific  heats  of  the 
species  (ftVsec2  6R) 

(INPUT) 

XW 

is  an  array  of  N molecular  weights  of 
the  species  (slug/slug  mole) 

(INPUT) 

T 

is  the  temperature  of  the  gas  fR) 

(INPUT) 

C 

is  an  array  of  N mass  fractions  of  the 
species 

(INPUT) 

IT 

is  an  array  of  N indices  of  the  species 
into  a table  of  collision  diameter  (cr) 
and  energy  of  attraction  (e/k) 

(INPUT) 

N 

is  the  number  of  species  in  the  gas 

(INPUT) 

XM 

is  an  array  of  N viscosities  of  the 
species,  (Ibf • sec/ft2) 

(0UTPUT) 

XMT 

is  the  total  viscosity  of  the  gas, 
(lbf  • sec/ft2) 

(0UTPUT) 

TK 

is  the  thermal  conductivity  of  the  gas, 
(ft.lbf/ft2  sec  PR/ft)) 

(0UTPUT) 

PR 

is  the  Prandtl  number  of  the  gas, 

(0UTPUT) 

Method: 

N 

C ^ E C.  C Specific  Heat 

p i=i  1 pi 


M 


w 


-fir 

2 

i=i 


Molecular  Weight 
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M 


» C.  . 


w 


i * M 


w, 


*1 


T/1.8 

urn; 


fl  = table  (T*) 


table  of  n vs  T* 


= table  (i) 


«\/M  T 

4.15822x10"  V wi 

/ M \ — 1/2 
1 ll+  Wl 


n r / n x.  \ 
= Si  Ui  (1+5i  xij 

L X tfi  / 


table  of  a and  c/k 
vs,  individual 
species 


».  \ ^ i \ ' 1/1 


J 


-1 


W 


W, 


viscosity  of  the  gas 


K, 


K 


M,  R 


M 


(.45  + 1.32 


Pi 


wx 


rnr) 


) 


/ N 

!■  K£  f 1 + lo065  S 


N 

r. 

1=1 


j.) 

i i xi 1 


-i 


j=i  i 

i^i 


thermal  conductivity 


K 


Prandtl  number 
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Equations  for  and  ]Ui  are  from  Reference  7.  The  values  of  the  collision 
integral  are  from  table  2 of  Appendix  B in  Reference  8.  The  relations  used 
to  calculate  /land  K of  the  mixture  are  from  Reference  9 . 

'Also  from  Reference  7 are  the  values  of  the  collision  diameters,  a,  and  energy 
of  attraction,  e/k. 

Table  1 correlates  the  chemical  name  of  the  species  to  the  internal  number 
assigned  to  it  by  the  subroutine.  Also  included  in  Table  1 is  the  key-punch 
name  assigned  to  'die  species,  since  lower-case  letters  and  subscripts  are 
non-standard  features  in  most  computer  configurations. 
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Table  1 . SPECIES  NAMES  AND  IDENTIFIERS 


Species 

Number 

Chemical 

Name 

Key  Punch 
Name 

1 

A1 

AL 

2 

A1C1 

ALCL 

3 

AlCla 

ALCL3 

4 

A1F 

ALF 

5 

A1F3 

ALF3 

6 

AIN 

ALN 

•/ 

A10 

ALO 

3 

A1S 

ALS 

9 

Ala 

AL2 

10 

Air 

AIR 

11 

Ar 

AR 

12 

AsH3 

ASH3 

13 

B 

B 

14 

BBr3 

BBR3 

15 

BC1 

BCL 

10 

BC12 

BCL2 

17 

BC13 

BCL  3 

16 

BF 

BF 

19 

bf2 

BF2 

20 

bf3 

BF3 

21 

bi3 

BI3 

22 

BO 

BO 

23 

b(och3)3 

B(OCH3)3 

Species 

Number 

Chemical 

Name 

Key  Punch 
Name 

35 

Br 

BR 

36 

BrF 

BRF 

37 

BrF3 

BRF3 

38 

BrO 

BRO 

39 

Br2 

BR2 

40 

C 

C 

41 

CBrF3 

CBRF3 

42 

CBr4 

CBR4 

43 

CC1 

CCL 

44 

CC1F3 

CCLF3 

45 

CC12 

CCL2 

46 

CC12F2 

CCL2F2 

47 

CC13 

CCL3 

48 

CC13F 

CCL3F 

49 

ecu 

CCL4 

40 

CF 

CF 

51 

cf2 

CF2 

52 

CF:j 

CF3 

53 

cf4 

CF4 

54 

CH 

CH 

55 

CHBrCIF 

CHBRCIF 

56 

CHBrCU 

CHBRCL2 

57 

CHBr3 

CHBR3 

58 

CHC1F2 

CHCLF2 

59 

CHClg 

CHCL3 

GO 

chf3 

CHF3 

Cl 

CHpBrCl 

CH2BRCL 

62 

CHnCIF 

CH2CLF 

63 

CHSCL«, 

CH2CL2 

64 

CH2F2 

CH2F2 

G5 

ch3i2 

CH2I2 

66 

CH3Br 

CH3BR 

67 

CH3CI 

CH3CL 
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Chemical 

Name 


Species 

Number 

58 

69 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79 

80 
81 
82 
03 

84 

85 

86 
37 
88 

89 

90 

91 

92 

93 

94 

95 
36 

97 

98 

99 


CH3F 

CH3I 

CH3OH 

CH4 

CN 

CO 

cos 

co2 

CP 

cs 

CSa 

Cs 

QjHa 

CSH 

CsHq 

CjHgCl 

CsHgOH 

C3N2 

CH3OCH3 

ch3chch3 

CH3CCH 

cyclo-CsHg 

C3He 

n-CsH^OH 

CH3COCH3 

CH3COOCH3 

n-C4H10 

iso-C4H10 

CsHsOCbHs 

CHsCOOCsHg 

n-C5H13 

C(CH3)4 


Key  Punch 
Name 

CH3F 

CH3I 

CH30H 

CH4 

CN 

CO 

COS 

C02 

CP 

CS 

CS2 

C2 

C2H2 

C2H4 

C2H6 

C2H5CL 

C2H50H 

C2N2 

CH30CH3 

CH2CHCH3 

CH3CCH 

CYCLO-C3H6 

C3H8 

N-C3H70H 

CH3COCH3 

CK3C00CH3 

N-C4H10 

ISO-C4H10 

C2H50C2H5 

CH3COOC2H5 

N-C5H12 

C(CH3)4 


1 


Species 

Chemical 

Key  Punch 

N amber 

Name 

Name 

101 

C6H12 

C6H12 

102 

n-CeHx* 

N-C6H14 

103 

Cd 

CD 

104 

Cl 

CL 

105 

C1CN 

CLCN 

106 

C1F 

CLF 

107 

C1F3 

CLF3 

108 

CIO 

CLO 

109 

Cl2 

CL  2 

110 

F 

F 

111 

FCN 

FCN 

112 

f2 

F2 

113 

H 

H 

114 

HBr 

HBR 

115 

HCN 

HCN 

116 

HC1 

HCL 

117 

HF 

HF 

118 

HI 

HI 

119 

HS 

HS 

120 

h2 

K2 

121 

h2o 

H2G 

122 

h2o2 

H202 

123 

1 .-sb 

H23 

124 

He 

HE 

125 

Hg 

HG 

126 

Hg  Br2 

HGBR2 

127 

HgCle 

HGCL2 

128 

Ugh 

HGI2 

129 

I 

I 

130 

ICI 

ICL 

131 

h 

12 

132 

Kr 

KR 

133 

Li 

LI 

134 

LiBr 

LIBR 
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Species 

Chemical 

Key  Punch 

Number 

Name 

Name 

135 

LiCN 

LICN 

136 

LiH 

LI  CL 

137 

LiF 

LIF 

138 

Lil 

LII 

139 

LiO 

LIO 

140 

\j1l2 

LI2 

141 

LijjO 

1/20 

142 

Mg 

MG 

143 

MgCl 

MGCL 

144 

MgCla 

MGCL2 

145 

MgF 

MGF 

146 

MgF2 

MGF2 

147 

Mg2 

MG2 

148 

N 

N 

149 

NFa 

NF3 

150 

NH 

NH 

151 

nh3 

NH3 

152 

NO 

NO 

153 

NOC1 

NOCL 

154 

N 

2 

N2 

155 

n2o 

N20 

156 

Na 

NA 

157 

NaBr 

NABR 

158 

NaCN 

NACN 

159 

NciJl 

NACL 

160 

NaF 

NAF 

161 

Nal 

NAI 

162 

NaO 

NAO 

163 

NaOH 

NAOH 

164 

Naa 

NA2 

165 

Na20 

NA20 

166 

Ne 

NE 

167 

0 

0 

168 

OF 

OF 
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Species 

Chemical 

Key  Punch 

Number 

Name 

Name 

169 

OFa 

OF2 

170 

OH 

OH 

171 

Oa 

02 

172 

P 

P 

173 

PCI 

PCL 

174 

PC13 

PCL3 

175 

PF 

PF 

176 

pf3 

PF3 

177 

ph3 

PH3 

178 

PN 

PN 

179 

PO 

PO 

180 

PS 

PS 

181 

P2 

P2 

182 

P4 

P4 

183 

S 

S 

184 

sf3 

SF6 

185 

SO 

SO 

186 

so2 

S02 

187 

S2 

S2 

188 

S2F2 

S2F2 

189 

SI 

SI 

190 

SiCl 

SICL 

191 

SiCk 

SICL  4 

192 

SlF 

SIF 

193 

SiFCla 

SIFCL3 

194 

SiF Cl 

2 2 

SIF2CL2 

195 

SiF.Cl 

SIF3CL 

196 

SiF4 

SIF4 

197 

SIH*. 

SIH4 

198 

SiO 

SIO 

199 

Si02 

SI02 

200 

SiS 

SIS 

°01 

2 

SI2 

202 

SnBr2 

SNBR2 

SnCU 

SMCL4 

204 

UFt, 

UF6 

205 

Xo 

XE 

206 

Zn 

ZN 
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5.3.12  Subroutine  C'DES 

This  is  the  m*tn  program  for  0DE  and  corresponds  to  the  Main  Program  de- 
scribed in  Reference  1.  Generally,  the  routine  oerforms  the  following  functions: 

1.  Reads  code  cards  THERM0,  REACTANTS,  0MIT,  INSERT,  and  NAMELISTS 
and  directs  flow  of  program  accordingly. 

2.  Stores  THERM0  data  on  tape. 

3.  Calls  subroutine  REACT  to  read  and  process  REACTANTS  cards. 

4.  Reads  0MIT  and  INSERT  cards  and  stores  species  names. 

5.  Initializes  variables  in  namelist  $0DE. 

6.  Reads  and  writes  namelist  $0DE. 

7.  Converts  assigned  densities,  if  any,  (RHO(i)  in  $0DE)  to  specific 
volumes:  VLM(i)  = l/RHO(i). 

8.  Stores  the  number  of  pressures  or  volumes  in  NP. 


9.  Stores  values  of  o/f  in  0XF  array.  If  o/f  values  have  not  been  input 
directly,  they  are  calculated  as  follows: 


Values 

Code 

o/f  calculation  in 
main  program 

Oxidant  to  fuel  weicnt  ratio,  o/f 
Fuel  to  air  weight  ratio,  f/a 
Percent  fuel,  %F 

Equivalence  ratio,  r 
Not  specified 

OF  = .TRUE. 

FA  = .TRUE. 
FPCT  = .TRUE. 

ERATI0  = .TRUE. 

0XF(i)  = o/f 
0XF(l)  = 1/  (f/ a) 

0XF(i)  = (100-%F)/(%F) 

_rV-(2)  v+(2) 
rtvp/i)  _ Tv  “V 

V^v1/  _(1)  + (1) 

rV  U'+V  u; 

rtvrfO  - WPW 
0XF(l)  WP(2) 

Values  of  WP(1)  and  V/P(2)  are  defined  in  appendix  B or  Reference  1. 

10.  Makes  necessary  adjustments  to  consider  charge  balance  if  I0NS= 
.TRUE. . This  is  done  by  adding  1 to  NLM  and  E to  LLMT  array. 

11.  Calls  SEARCH  to  pull  required  THERM0  data  from  tape  and  to  store  the 
data  in  core. 

12.  Sets  initial  estimates  for  compositions.  These  estimates  are  set  with 
each  $0DE  read.  They  are  used  only  for  the  first  point  in  the  lists  of 
variables  in  namelist  (e.g.,  the  first  o/f  and  the  first  T and  P in  a TP 
problem).  /II  succeeding  points  use  results  from  a previous  point  for 
estimates. 
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For  the  first  ooint  the  program  assigns  an  estimate  of  0.1  for  n,  the  total 
number  of  kilogram  -moles  per  kilogram.  The  initial  estimate  of  number  of  moles 
of  each  gaseous  species  per  kilogram  of  mixture  n^  is  set  equal  to  0.1/m  where 
m is  the  total  number  of  gaseous  species.  Condensed  species  are  assigned  zero 
moles . 

13.  Sets  IUSEft)  positive  for  condensed  species  listed  on  INSERT  cards 
(see  IUSF.  array). 

14.  Calls  R0CKET  is  RKT  is  true. 

Note:  while  the  use  of  oxidizer  to  fuel  ratio,  0/F,  is  not  generally  used 
in  describing  solid  rocket  motor  propellants,  it  can  be  used  to  great  advantage  in 
specifying  binder,  metal  loading,  and  oxidizer  ratios  in  parametric  studies. 


5.3,13  Subroutine  OMEGA 
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This  subroutine  calculates  the  exponent,  used  in  the  viscosity-temper- 
ature relationship 

u=Mref  ( 


T 60 

r~  ) 


ref 


using  the  method  of  least  squares.  That  is,  it  calculates  the  value  of  fJt3  which 
gives  the  smallest  sum  of  the  errors  squared.  This  form  of  the  viscosity- temper- 
ature relationship  was  selected  since  both  the  TD2P  and  TBL  modules  require 
viscosity  data  in  this  manner. 

In  order  to  supply  the  maximum  amount  of  accuracy  and  also  to  minimize  the 
variation  in  data  due  to  the  selection  of  an  exit  area  ratio,  it  was  decided  to  match 
the  throat  value  of  viscosity  exactly  and  select  an  which  would  provide  the  best 
fit  for  viscosity  at  the  chamber  and  exit  of  the  motor. 

The  form  of  the  error,  E,  was  taken  to  be 


E = In  a/fi*  " co  ln  T/T* 

Squaring  the  errors,  differentiating  with  respect  to  oo  , and  setting  the  results 
equal  to  zero,  yields  the  following  value  for  ^ 

os  = (In  Tc/T*  In  4/M*  + In  T/T*  In  nf/fl*)Aln  T/T*)2  + (In  T/T*)2) 


where 


T = temperature 
H = viscosity 
e = refeis  to  the 
c .=  refers  to  the 
* = refers  to  the 


exit  plane 
chamber 
throat  plane 
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5.3.14  Subroutine  0UT1 

This  subroutine,  together  with  entries  0UT2  and  0UT3,  writes  statements 
common  to  all  problems.  0UT1  writes  statements  giving  the  data  on  REACTANTS 
and  on  o/f,  percent  fuel,  equivalence  ratio,  and  density. 

Entry  0UT2 . - This  entry  writes  the  statements  for  printing  values  of  pres- 
sure, temperature,  density,  enthalpy,  entropy,  molecular  weight,  fe  In  V/d  In  P),^ 
(if  equilibrium),  fe  In  V/fe  In  T)p  (if  equilibrium),  heat  capacity,  yg,  and  sonic 
velocity.  These  variables  and  corresponding  labels  are  printed  with  a variable 
format  described  in  BL0CK  DATA. 

Entry  0UT3.  - Entry  0UT3  writes  statements  giving  the  equilibrium  mole 
fractions  of  reaction  species. 

Subroutine  0 CTT1  has  also  been  modified  to  print  out  the  following  addition- 
al items: 

a)  gas  velocity 

b)  molecular  weight  of  the  gas-condensed  phase  combined 

c)  mass  fraction  of  condensibles 

d)  molecular  weight  of  the  condensibles  in  the  chamber 

e)  molecular  weight  of  the  gas  in  the  chamber 

f)  the  erosion  parameter, 

In  addition,  0UT1  also  traps  the  starting  conditions  for  the  0DK  module. 
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5.3.15  Subroutine  REACT 


The  purpose  of  subroutine  REACT  Is  to  read  and  process  the  data  on  the 
REACTANTS  cards.  The  subroutine  is  called  from  the  main  program  after  a REACTANTS 
code  card  has  been  read.  The  data  on  these  cards  are  described  in  the  REACTANTS 
Cards  section  (p.  62)  of  Reference  1 and  in  Volume  III,  Section  2.6.1.  References 
to  page  numbers  and  equations  given  below  also  pertain  to  Reference  1 . 

The  reactants  may  be  divided  into  two  groups  according  to  card  column  72 
on  th?  REACTANTS  cards.  The  two  groups  are  oxidants  (O  in  cc  72)  and  fuels  (cc 
72  £ O).  We  generally  keypunch  F in  card  column  72  for  fuels  even  though  this 
is  not  necessary.  The  contents  of  card  column  72  are  read  into  F0X.  Depending 
on  the  contents  of  F0X,  program  variables  relating  to  oxidants  or  fuels  are  sub- 
scripted 1 for  oxidants  and  2 for  fuels. 

The  FORTRAN  symbols  for  the  properties  read  from  the  REACTANTS  cards 
and  their  associated  properties  (discussed  in  INPUT  CALCULATIONS,  p.  55  of  Ref- 
erence 1)  are  as  follows: 


Property 

FORTRAN  symbol 

a(k) 

aij 

wjk) 

Njk) 

(H°  )(k) 

(u^)(k) 

j 

00 

pj 

ANLM(j,m)a 

PECWT(j)  (if  no  M in  cc  53) 

PECWT(j)  (if  M in  cc  53) 

ENTH(j)  (if  not  UV  Problem  and  00  not  in  cc  37  and  38) 
ENTH(j)  (if  UV  problem  and  00  not  in  cc  37  and  38) 

DF,NS(j) 

1 

aEach  of  the  j REACTANTS  cards  contains  from  1 to  5 stoichiometric  coefficients 
read  (indicated  by  subscript  m)  into  ANUM(j,m)  and  their  corresponding  chemical 
symbols  read  into  NAME(j,m).  In  relating  an  ANUM(j,m)  with  afj^  , the  index  i 
associated  with  a particular  chemical  element  is  determined  from  the  chemical 
symbol  in  NAME(j,m). 
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If  there  are  several  oxidants  their  properties  are  combined  by  subroutine 
REACT  into  properties  of  a total  oxidant  using  the  relative  proportion  of  each 
oxidant  given  on  the  REACTANTS  cards*  Similarly,  if  there  are  several  fuels, 
their  properties  are  combined  into  properties  of  a total  fuel.  The  total  oxidant 
and  total  fuel  properties  discussed  in  INPUT  CALCULATIONS^  and  their  associated 
FORTRAN  symbols  are  as  follows: 


Property 

FORTRAN  symbol 

Equation 
(see  Ref,  1) 

b[k> 

BCP(i,k) 

(187) 

Rr<fW(j) 

(190) 

v(k) 

HPP(k)  (if  not  UV  problem  and  00  not  in  cc  37  and  38) 

(192) 

HPP(k)  (if  UV  problem  and  00  not  in  cc  37  and  38) 

(194) 

AM(k) 

(196) 

(k) 

0 

RH(k) 

(198) 

V+(  k) 

VPLS(k) 

(200) 

v”  M 

VMIN(k) 

(201) 

PECWT(j) 

If  any  of  the  p are  zero  then  RH(1)  = RH(2)  = 0. 

These  total  oxidant  and  total  fuel  properties  are  subsequently  combined 
into  total  reactant  properties  by  using  the  values  of  oxidant-fuel  mixture  ratios 
obtained  from  the  main  program.  This  Is  done  in  NEW0F,  an  entry  point  in  SAVE. 

Other  common  variables  set  by  REACT  are  LLMT,  NAME,  ANUM,  ENTH,  FAZ, 
RTEMP,  F0X,  DENS,  RMW,  M0LES,  NLM,  NEWR,  and  NREAC. 

A provision  is  made  for  eliminating  a second  tape  search  when  two  consecu- 
tive sets  of  REACTANTS  cards  contain  the  same  elements.  This  is  done  by  saving 
the  element  symbols  (LLMT(  ?))  in  LLMTS(i),  the  kilogram-atoms  per  kilogram 
(BOPU  ,k))  in  SB0PU,k),  and  the  number  of  elements  (NLM)  in  NLS, 

Atomic  weights  used  in  equation  (190)^  are  stored  in  AT0M(2,i).  The 
corresponding  chemical  symbols  are  stored  in  AT0M(l,i).  The  oxidation  states 
of  the  chemical  elemenis  V*or  V~  used  in  equations  (200)  and  (201)^  are  stored 
in  AT0M(3,i),  The  AT0M  array  is  stored  by  BL0CK  DATA. 
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5,3. 16  Subroutine  RKT0UT 


This  subroutine  calculates  various  rocket  performance  parameters  from  pre- 
viously calculated  thermodynamic  properties. 

It  is  also  the  control  program  for  writing  rocket  performance  output.  It  con- 
tains the  WRITE  statements  that  apply  specifically  to  rocket  parameters  and  it 
calls  subroutine  0TJT1  and  entries  0UT2  and  0UT3  for  the  WRITE  statements  common 
to  all  problems.  The  rocket  parameters  are  printed  with  the  variable  format,  FMT, 
described  in  BL0CK  DATA. 

Subroutine  RKT0UT  is  called  from  subroutine  R0CKET. 
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5.3.17  Subroutine  R0CKET 

This  subroutine  is  the  control  program  for  the  RKT  problem  (rocket  perfor- 
mance calculations  discussed  in  section  R0CKET  PERFORMANCE).^  A flow 
diagram  for  this  subroutine  is  given  in  Figure  5 of  Reference  10  Subroutine 
R0CKET  obtains  the  required  thermodynamic  properties  for  equilibrium  performance 
by  calling  subroutine  EQLBRM.  For  forzen  performance,  subroutine  R0CKET  calls 
subroutine  FR0ZEN  to  obtain  the  required  thermodynamic  properties.  Rocket  per- 
formance parameters  are  then  obtained  by  calling  subroutine  RKT0UT.  In  addition 
to  calling  RKT0UT,  and  FR0ZEN,  and  in  addition  to  using  controls  common  to  all 
problems  (discussed  in  section  MODULAR  FORM  OF  THE  PROGRAM,  p.  75,  Ref- 
erence 1)  subroutine  R0CKET  also  does  the  following: 

1.  It  calculates  estimates  for  throat  pressure  ratios. 

2.  It  calculates  estimates  for  pressure  ratios  corresponding  to  assigned 
area  ratios  (if  any). 

Subroutine  R0CKET  was  modified  to  perform  the  following  tasks: 

a)  compute  the  area  ratio  at  which  solidification  of  condensed  phases 
would  occur  if  condensed  phases  are  present. 

b)  on  option,  restricts  the  total  amount  of  condensed  phases  during 
the  expansion  to  the  amount  Initially  in  the  rocket  motor  chamber. 

c)  computes  the  print  out  station  at  which  subroutine  OUT1  will  trap 
the  start  conditions  for  a kinetics,  0DK,  calculation. 

d)  computes,  using  subroutiens  MATCH,  MUK,  and  0MEGA,  transport 
properties  for  the  BAL,  TD2P,  and  TBL  modules. 

e)  writes  (punches)  out  the  linkage  data  necessary  for  communications 
with  the  various  loss  modules. 
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5,3,18  Subroutine  SAVE 

This  subroutine  has  several  functions,  all  of  which  are  concerned  with  sav- 
ing some  information  from  a completed  calculation  for  subsequent  use  in  later  cal- 
culations, The  primary  purpose  is  to  save  computer  time  by  having  good  initial 
estimates  for  compositions. 

These  estimates  for  the  next  point,  NPT,  come  from  either  the  point  just 
completed,  ISV,  or  some  other  previous  point.  The  flow  of  the  routine  is  directed 
by  ISV  as  follows: 

1.  ISV  positive.  Transfer  compositions  from  the  point  just  completed  for 
use  as  initial  estimates  for  next  point  (transfer  EN(j,  ISV)  to  EN(j,  NPT)), 

2.  ISV  negative.  Save  values  of  ENLN(j)  for  gases  and  EN(j)  for  condensed 
in  SLN(j) , ENNin  ENSAVE,  ENNL  in  ENLSAV,  IQ1  in  IQSAVE,  JS0L  in 
JS0LS,  JL1Q  in  LKIQS,  and  NLM  in  LL1.  (These  values  are  saved  because 
they  are  to  be  used  as  Initial  estimates  for  some  future  point  and  they  may 
be  overwritten  in  the  meantime.)  Make  ISV  positive  and  transfer  EN(J,ISV) 
to  EN(j,NPT). 

3.  ISV  zero.  Use  the  data  previously  saved  (as  discussed  in  2.  as  initial 
estimates  for  current  point.  Restore  IUSE  codes  and  inclusion  or  exclu- 
sion of  "E"  as  an  element  for  I0NS  option. 

Entry  NEW0F.  - NEW0F  combines  the  properties  of  total  oxidant  and  total 
fuel  calculated  in  subroutine  REACT  with  an  o/f  value  to  give  properties  for  the 
total  reactant.  NEW0F  is  called  for  each  mixture  assigned  in  the  MIX  array  in 
$0DE  namelist.  It  is  called  from  subroutine  R0CKET.  The  calculated  properties 
and  corresponding  FORTRAN  symbols  are  as  follows: 


Property 

FORTRAN  symbol 

b°i 

h(/R 

u^/R 

po 

r 

B0(l) 

HSUBO  (If  not  UV  problem) 
HSUBO  (If  UV  problem) 
RH0P 
EQRAT 

Equation 


(191) 

(193) 

(195) 

(199) 

(204) 


Subroutine  HCALC  Is  called  by  Entry  NEW0F  to  calculate  the  enthalpies 
for  each  reactant  that  has  zeros  keypunched  In  card  columns  37  and  38  In  its 
REACTANTS  card. 


Values  of  HPP(2),  HPP(l).  HSUBO,  B0P(l,2),  B0P(t,l),  and  B0(l)  are  printed 
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5.3.19  Subroutine  SEARCH 


This  subroutine  selects  the  Thermodynamic  Data  to  be  used  in  the  problem. 

A scan  is  made  of  the  master  Thermodynamic  Data  tape  and  those  species  that  are 
consistent  with  the  chemical  system  under  consideration  are  selected.  As  the 
thermodynamic  data  are  being  selected,  the  subroutine  also  complies  a set  of 
formula  numbers,  a^,  from  the  formulas  of  the  reaction  products.  A short  Thermo- 
dynamic Data  file  is  also  generated  for  use  in  subsequent  calculations, 

A check  is  made  near  the  beginning  of  the  routine  to  prevent  THERM0  data 
from  exceeding  their  storage  allotments.  These  variables  are  all  in  labeled  common 
SPECIES  and  are  currently  dimensioned  for  150  species.  However,  this  dimension 
may  be  reduced  to  save  storage. 

SEARCH  is  called  from  the  main  program  when  the  logical  variable  NEWR 
is  true.  NEWR  is  set  true  in  REACT  to  indicate  a new  chemical  system.  REACT 
also  stores  chemical  element  symbols  for  the  current  chemical  system  in  the  LLMT 
array.  SEARCH  stores  THERM  0 data  in  core  for  each  species  whose  elements  are 
included  in  the  LLMT  array  (unless  the  species  name  was  listed  on  an  0MIT  card). 

The  THERM0  data  are  stored  in  common  variables  TL0W,  TMID,  THIGH, 

SUB,  A,  C0EF,  and  TEMP.  SEARCH  writes  out  the  names  and  dates  of  species 
whose  data  are  stored  in  core. 

SEARCH  initializes  the  IUSE  array.  IUSE(j)  for  gaseous  species  are  set 
equal  to  zero.  IU3E(j)  for  condensed  species  are  set  equal  to  negative  integers. 

For  the  chemical  system  under  consideration  the  first  possible  condensed  species 
is  set  equal  to  -1,  the  second  to  -2,  and  so  on,  with  one  exception.  In  the  event 
there  are  two  or  more  condensed  phases  of  the  same  species,  set  equal  to  -4,  for 
example,  IUSE(j)  for  B^O^s)  will  also  be  set  equal  to  -4.  A description  of  the  IUSE 
array  is  given  below. 

The  various  condensed  phases  of  a species  are  expected  to  be  adjacent  in 
the  THERM 0 data  as  they  are  read  from  tape.  These  phases  must  be  either  in 
increasing  or  decreasing  order  according  to  their  temperature  intervals. 

NS  contains  the  total  number  of  species  stored  in  core.  NC  contains  the 
total  number  of  condensed  species  (counting  each  condensed  phase  of  a species 
as  a separate  species). 
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IUSE  array.  - Each  value  in  the  IUSE  array  Is  associated  with  a species. 
These  values  of  IUSE  serve  two  purposes: 

1.  They  indicate  which  species  are  to  be  included  in  the  current  iteration 
(IUSE (j)  <0  for  excluded  species  and  IUSE(j)  ^0  for  included  species). 

2.  They  indicate  multiple  phases  of  the  same  species  if  absolute  values  of 
IUSE(j)  are  equal. 

The  IUSE(j)  are  initialized  in  subroutine  SEARCH  and  the  main  program  as 
follows: 

1.  IUSE(J)  = 0 for  all  gaseous  species. 

2.  IUSE(j)  = n for  all  condensed  species  whose  names  have  been  listed  on 
INSERT  cards.  The  number  n indicates  the  species  was  the  n*h  condensed 
species  whose  THERM  0 data  were  read  from  tape. 

3.  IUSE(j)  = ~n  for  all  condensed  species  not  listed  on  INSERT  cards  where 
n is  defined  in  2. 

These  initial  values  of  IUSE(j)  may  be  adjusted  later  in  subroutine  EQLBRM. 
For  condensed  species,  the  sign  is  adjusted  as  species  are  included  or  excluded 
in  the  current  iteration. 

For  the  I0NS  option,  IUSE(j)  values  for  ionic  species  are  set  to  -10000  when 

_o 

the  mole  fractions  of  all  ionic  species  are  less  chan  10 


5.3,20  Subroutine  SHCK 

Subroutine  SHCK  Is  a dummy  routine  which  replaces  subroutine  SHCK  of 
Reference  1. 


5.3.21  Subroutine  THERM P 

Subroutine  THERMP  Is  a dummy  routine  which  replaces  subroutine  THERM P 
of  Reference  1. 


5.3.22  Subroutine  VARFMT 

Subroutine  VARFMT  (variable  format)  adjusts  the  number  of  decimal  places 
printed  in  F-format  in  the  variable  format,  FMT,  according  to  the  size  of  the 
number.  It  Is  used  for  P </?  , P,  and  Ae/At-  Variable  format  is  described  in 
BL0CK  DATA , Section  5.1.2. 


5,4  Link  30  Subroutines 


The  Link  30  subroutines  control  the  grain  design  (GD)  and  ballistics  (BAL) 
calculations  made  in  the  SPP  code.  These  two  modules  have  been  integrated  to- 
gether from  the  programs  of  Reference  2 and  Reference  6.  The  remainder  of  this 
subsection  gives  an  overall  description  of  the  integrated  ballistics  module. 
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5.4,0, 1 General 


The  Interior  Ballistics  Module  models  the  Internal  ballistics  of  solid  rocket 
motors.  Including  motor  environmental  effects  on  burning  rate  (erosion,  radiation), 
effects  of  radial  slots  between  segments  on  the  port  flow,  port  flow  pressure  drop, 
and  delayed  burning  (ignition)  in  deep  narrow  slots.  Semi-empirical  correlations 
for  combustion  inefficiency  and  nozzle  throat  erosion  also  are  included.  This 
computer  module  was  adapted  from  an  existing  model  used  to  analyze 
nozzleless  solid  rocket  motors  (Ref.  6),  because  it  has  gasdynamic  features 
lacking  in  many  published  internal  ballistics  programs.  The  intent  was  to  predict 
internal  ballistics  for  motors  ranging  from  small  tactical  to  large  boosters  wherein 
aspects  of  these  features  are  encountered.  The  major  effort  in  the  present  program 
development  Involved  updating  certain  model  expressions,  and  rendering  the  com- 
putational procedure  compatible  with  the  Grain  Design  Module  and  the  rest  of  the 
SPP  code.  The  model  will  be  subject  to  continued  revision  as  new  information  be- 
comes available. 

The  logic  of  solving  the  ballistics  equations  is  presented  in  the  following 
subsection.  A discussion  of  the  module  subroutine  structure  is  presented  after  that. 
Analogous  discussion  of  the  grain  design  subroutine  structure  is  contained  in  Ref. 

(2)  and  need  not  be  repeated  here.  Finally,  the  details  of  program  Input  and  out- 
put, for  the  integrated  subprogram,  are  given  in  Volume  III,  The  Program  User's 
Manual. 

The  overall  geometry-aerodynamics  interface  between  the  GD&BM  is  the 
generation  of  port  area  change  as  a function  of  length  in  the  Grain  Design  routines 
and  return  of  web  burned  from  BM  to  GD  module,  at  the  same  length  stations. 
Ballistic  calculations  may  be  performed  for  one  or  more  time  increments  within  an 
increment  of  bum  (NB)  specified  in  the  GD  input,  but  as  the  increment  is  predicted 
to  bum  out  the  actual  web  expended  is  returned  as  the  starting  point  for  the  next 
increment  of  geometry  calculation. 
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5«4»0.2  Solution  Logic 


The  general  scheme  employed  in  solving  the  equations  of  one-dimensional 
flow  is  shown  in  simplified  terms  in  Figure  5.4-1.  The  procedure  that  follows 
the  initialization  of  a case  is  as  follows: 

• A head-end  pressure  is  assumed.  The  calculations  for  a single  axial  in- 
crement are  made  and  iterated  to  completion.  The  successive  increments 
down  the  port  are  then  calculated  until  a match  point  on  the  nozzle  entrance 
is  reached.  * 

• At  the  match  point,  the  local  static  pressure  is  checked  against  the  pres- 
sure from  throat  flow  constraints.  If  the  two  values  are  not  within  tolerance, 
a correction  to  head-end  pressure  is  made  and  the  calculation  is  restarted 

at  the  head.  As  successive  iterations  are  made  to  satisfy  the  choke  con- 
straint (outer  loop),  the  axial  stations  (inner  loop)  are  continuously  iterated. 

• Parameters  used  in  transient  terms  and  integrating  calculations  are  saved, 
and  the  time  increment  is  advanced.  The  calculations  resume  with  a new 
head-end  pressure  assumption. 

• As  the  calculations  reach  the  end  of  a specified  run  time,  or  the  motor 
has  burned  out,  output  is  furnished  to  the  Nozzle  Performance  Module. 

The  computational  logic,  thus,  is  arranged  as  a triple-loop  structure  with 
the  axial  increment  (inner  loop),  choke  constraint  (outer  loop)  and  time  increments 
(stepping). 

A more  detailed  diagram  of  program  logic  is  shown  in  Figure  5.4-2.  This 
representation  is  still  simplified  with  regard  to  the  individual  branching  for  special 
cases,  but  is  intended  to  describe  the  function  of  groupings  of  program  steps  in- 
ternal to  the  program,  and  presumes  valid  input.  Figure  5.4-2  is  organized  in  the 
same  order  as  the  source  program. 

Ihe  compiled-in  constants  are  established  and  some  core  areas  are  cleared 


to  zero. 


A subroutine  to  control  the  reading  and  organizing  of  tabular  data  is  called 
b\  Subroutine  READIN  on  completion  of  geometry  inputs.  The  remaining  data  (in- 
tegers and  real  data  ir.out)  are  returned  in  intermediate  storage  which  is  inspected 
for  non-zero  (nz)  values.  Any  nz  values  are  then  placed  in  working  storage  for 
use  hi  MAIN  3. 


For  an  end-burner,  the  ‘head  end"  will  be  the  grain  face.  A protective  statement 
was  included  providing  that  if  the  motor  port/throat  is  less  than  1/4,  the  calcula- 
tion will  proceed  to  the  next  increment.  Thus,  the  end  face  is  eventually  located. 
Internal  burners  are  unlikely  to  have  such  low  port/throat,  so  this  provision  should 
not  introduce  a problem  in  internal  burners. 
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Hie  outside  diameter  or  periferal  geometry  of  the  grain  is  established  from 
tabular  inputs  which  are  fed  to  subroutine  L00PS  in  the  GD  program.  The  GD  pro- 
gram initializes  the  grain  and  computer  bumback  geometires. 

The  following  groups  of  calculations  are  repeated  at  each  time  increment 
of  the  problem.  The  igniter  flow,  if  any  is  determined,  and  time  dependent  values 
of  Aj#  Ae,  thrust,  ambient  pressure  and  thrust  coefficient  are  determined  either 
from  tabular  input  or  internal  feedback.  Calculations  are  made  for  each  axial 
increment  in  the  motor.  The  axial  increments  correspond  to  the  GD  X-grid 
plus  an  aft  plenum,  plus  the  throat.  Trial  values  for  the  initial  time  step  are  used 
from  the  previous  axial  station;  all  other  trial  values  result  from  the  previous  iter- 
ation of  the  station.  Special  case  logic  for  ignition  options,  or  initial  time  are 
not  shown  in  the  figure  or  discussed  here.  (However,  they  should  be  recognized 
in  the  program  listing  as  IGN.GT.O,  or  HME.EQ.O.  branches.) 

The  port  geometry  is  calculated  by  the  GD  module  from  the  web  burned. 
Results  are  interpolated  within  the  bum  interval  to  determine  local  bum  volumes. 

• Base  strand  bum  rate  for  the  propellant  at  the  axial  station,  as  a function 
of  local  pressure,  is  found  from  an  input  table.  The  input  data  can  be 
keyed  to  use  log -log  interpolation.  The  motor  environment  contributions 
are  then  calculated.  7 7^  effects  are  Incorporated  by  changing  the  base  table. 

• If  appropriate,  logic  and  branching  related  to  radial  slots  is  performed. 

The  logic  path  returns  several  places  depending  on  the  interpretation  of 
radial  slot  flow. 

•Local  conditions  are  calculated  using  the  energy,  state,  and  continuity 
equations,  and  thermochemical  data. 

• Traps  are  executed  to  prevent  crossing  Mach  1,  and  the  momentum  equation 
is  used  to  solve  for  pressure.  The  set  of  equations  at  each  axial  station 

is  iterated  with  logic  for  subsonic  flow  or  Mach  1,  with  generation  of  an 
error  term  to  later  correct  head  pressure;  excess  iterations  will  be  shown 
as  output. 

• Looping  logic  allows  10  iterations  and  then  accepts  the  result  whether 
converged  or  not. 

• The  node  number  of  the  axial  increment  is  checked  to  branch  to  the  next 
node  or  to  throat  calculations. 

• Throat  constraints  are  calculated  from  mass  flow  rate,  nozzle  throat  area, 

C*,  etc. 

•Nozzle  entrance  static  pressure  is  calculated  from  throat  stagnation  pressure, 
assuming  an  isentropic  expansion,  to  compare  with  the  pressure  determined 
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by  Integration  down  the  port.  An  updated  head  pressure  Is  calculated  from 
the  error  In  nozzle  Inlet  static  pressure,  or  a term  Is  created  to  Indicate 
an  overchoked  prior  calculation. 

***  i^iroat  constraint  is  not  satisfied,  loops  are  counted  to  10.  From  10 
to  20  loops,  the  enror  Is  printed;  after  20  loops,  the  result  Is  accepted. 

The  printed  error  messages  Indicate  whether  the  calculations  are  oscillating 
or  are  overdamped  and  allow  the  user  to  determine  whether  the  run  Is  acceptable. 

• Values  used  In  transient  terms  or  Integrations  are  stored  for  use  In  the  next 
time  step.  Parameters  along  the  motor  are  printed.  The  burn-back  of  radial 
slots  and  wet  are  updated.  Axial  nodes  are  coded  out  Just  prior  to  bum- 
back  of  the  radial  slot  face  to  the  node  position.  For  radial  slot  combus- 
tion the  slot  ends  are  Input  to  the  GD  program  as  non-burning.  The 

bum-back  Is  then  explicitly  calculated  In  MAIN3  and  the  axial  Increment 
locations  are  modified. 

• Nozzle  exit  parameters  and  thrust  are  calculated  from  isentropic  expansion 
equations,  and  a series  of  running  Integrals  are  calculated  and  output. 

• The  values  of  time  and  Input  run  time  are  checked  for  completion.  The 
next  time  Increment  or  rewind  follows  this  branch. 

The  Indexing  scheme  employed  In  the  program  Is  shown  In  Table  5.4-1 
as  an  aid  In  following  the  details  of  the  logic.  This  logic  Is  designed  to  treat  non- 
perpendicular  bum-back  of  radial  slots;  but  Is  also  used  to  "key"  pressure  drop 

calculations  such  as  a sudden  expansion  or  non-dlffuslng  taper  (l.e.,  separated 
flow). 


Table  5.4-1  INDEX  SCHEME 


Symbol 


Description 


I 

Axial  node  location  of  calculation. 

I 

Prior  axial  node,  usually  J=I-1,  J>0,  but 
is  varied  to  skip  cast-off  nodes  (for  example 
J=I-4  if  three  nodes  are  burned  out). 

L 

Radial  slot  number,  1 at  head  end;  slot  re- 
lates to  following  segment. 

N 

Node  number  of  nczzle  entrance  (36  in 
sketch) . 

NN 

Throat  node. 

8 


Thrust  calculations  and  output 
Aero  calculations  and  output 
Ballistics  calculations  and  output 
Integrals  and  output 


Check  for  completion 
of  burn  Increment 

I 

Check  for  completion 
No  - Advance  time  Increment 
yes 

No 

yes 

RETURN  KBAL 


RETURN  TO  GD 


Figure  5.4-2  Module  Structure 
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Figure  5.4-3  Module  Logical  Structure 


Figure  5.4-4  Input  Groupings  Following  Grain  Geometry  Cards 


5,4 . 1  Program  LINK40 


This  subprogram  is  mainly  intended  to  facilitate  the  conversion  of  the  SPP 
code  overlay  structure  to  other  computer  systems.  Its  only  function  is  to  call  sub- 
routine BALIS. 

5.4.2  Subroutine  BALIS 

This  subroutine  reads  and  writes  the  namelist  data  set  $BAL.  It  then  calls 
subroutine  READIN  and  FIG1  to  read  in  and  then  calculate  the  grain  design  parameters. 
Control  of  the  ballistics  and  grain  design  calculations  sequences  are  then  handled 
in  subroutines  MAIN3,  PVS  and  L00PS.  MAIN3  performs  the  axial  looping  for  bal- 
listics calculations,  while  L00PS  does  the  same  for  the  grain  geometry  calculations. 


5.4.3  BL0CK  DATA 

This  routine  is  used  to  zero  out  the  common  blocks  used  in  this  module. 
The  KALPHA  array  in  labeled  common  C0DE  is  also  initialized. 

5.4.4  Subroutine  C0NTR1 

This  subroutine  has  been  written  for  use  in  the  Grain  Design  Program  as 
an  auxiliary  program  for  CONTR2.  The  routine  computes  the  initial  data  needed 
by  CONTR2  that  is  dependent  only  on  the  type  of  cone  (or  cylinder)  being  consid- 
ered and  independent  of  the  burn  distance,  X plane,  or  line  in  question, 

5.4.5  Subroutine  C0NTR2 

This  subroutine  has  been  written  for  use  in  the  Grain  Design  Program. 

The  purpose  of  the  routine  is  to  determine  the  points  of  intersection,  if  any, 
made  by  a line  and  a frustrum  of  a cone  (or  cylinder)  located  in  three  dimensional 
space. 

5.4.6  Subroutine  C3TEFF 

This  subroutine  calculates  the  empirical  c*  efficiency  described  in  Volume 
I,  Table  4-2. 
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5,4.7  Subroutine  DA  TAIN 


TWAIN  is  n random  rnnd  routine  tint  Interprets  input  of  selected  variables. 
UATA1N  operates  under  the  control  of  TABIN,  which  requests  the  reading  of  cards 
with  integers  or  real  data,  which  is  in  turn  called  from  Subroutine  READIN. 

5,4.8  Subroutine  ERROR 

This  subroutine  was  written  for  use  with  the  Grain  Design  Program.  The 
routine  checks  the  input  data  for  errors  and  writes  error  messages.  See  Appendix 
A for  conversion  aides. 


5.4.9  Subroutine  FIG1 

Subroutine  FIG1  has  been  written  for  use  with  the  Grain  Design  Program. 
The  program  calls  the  figure  routines  for  initial  calculations.  The  initial  bum 
used  to  input  geometric  figures  with  rounded  comers  is  monitored  by  subroutine 
FIG1. 


As  Subroutine  FIG1  is  entered,  I is  initialized  to  one  for  the  first  figure. 

ID  is  computed  from  IDENT(I)  and  a "computed  go  to"  is  executed  on  ID  in  order 
to  call  the  proper  figure  routine  (CONTRJ  or  PRISMl).  On  return  from  the 
figure  routine,  I is  incremented  by  the  proper  amount  and  a test  of  I against 
NOSRF  is  made  to  determine  if  more  figures  remain.  Upon  completion  of  the  figures, 
control  is  returned  to  the  main  program. 

5,4.10  Subroutine  FIG2 

Subroutine  FIG2  has  been  written  for  use  with  the  Grain  Design  Program. 

The  program  calls  the  proper  figure  routines  for  computing  the  points  defining  the 
void  volume  of  the  propellant. 

As  Subroutine  FIG2  is  entered,  the  motor  case  boundaries  are  computed  for 
the  specified  Y value  being  considered.  Then  the  bum  distance  to  be  considered 
(TRIAL)  and  the  identification  code  (ID)  for  the  figure  to  be  considered  are  com- 
puted. Then  the  proper  figure  routine  is  called  and  upon  return  Subroutine  UNIOr 
is  called  to  form  a union  between  the  computed  Z points.  These  points  represent 
the  shape  of  the  propellant  of  the  motor  case. 
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If  all  input  figures  have  not  beer  covered,  the  loop  is  repeated,  from  the  computing 
of  TRIAL,  for  the  next  input  figure. 

5.4.11  Subroutine  HNDX 

This  routine  performs  linear  interpolation  or  extrapolation  in  tables  while 
saving  its  place  in  the  table  of  independent  variables. 

5.4.12  Subroutine  L0(2flPS 

Subroutine  L00PS  sequences  the  grain  cfesicn  subroutines  through  the  input 
lists  of  Y increments  and  performs  a summation  cf  port  area  contributions  for  the 
Ntr  axial  station  in  the  motor.  This  subroutine  replaces  several  in  the  original 
program  of  Pef.  (2)  since  less  options  are  retained  for  ballistics  calculations  and 

inverted  order  of  sequencing  is  not  needed.  Thus,  in  the  revised  version,  the  outer 

? 

loop  is  always  bum  increment;  the  middle  loop  is  always  motor  axial  station  (di- 
mension x);  and  the  inner  loop  is  on  the  radial  intersection  (dimension  y).  The  sub- 
routines performing  the  actual  geometry  calculation  are  essentially  unchanged  from 
Ref.  (2):  changes  were  introduced  in  the  sequencing  routines  to  imporve  program 
run  times . 

5.4.13  Subroutine  MAIN3 

This  routine  controls  the  calculation  of  both  the  ballistics  and  grain  design 
modules.  In  fact  it  is  the  main  portion  of  the  ballistics  calculation. 

Subroutine  MAIN3  replaces  a series  of  subroutines  controlled  by  MAIN3  of 
the  original  Ref.  (2)  program.  This  subroutine  calculates  internal  motor  ballistics 
based  on  grain  geometry  output  from  the  grain  design  subroutines.  Control  is 
retained  in  MAIN3  fcr  as  many  time  increments  of  web  burned  as  necessary  to  pre- 
dict bumcut  of  the  particular  increment  of  web  calculated  by  the  grain  design  pro- 
gram. The  web  calculated  by  the  ballistics  subroutines  is  then  returned  to  the  grain 
design  program  through  subroutine  BALIS  for  use  as  the  starting  web  in  the  next 
Increment. 

An  artificial  point  in  the  nozzle  entrance  region  is  introduced  as  a match 
point  for  calculations  proceeding  down  the  port  to  the  nozzle,  for  the  purpose  of 
testing  the  choking  constraint  without  having  to  calculate  flow  near  Mach  one. 
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The  area  of  the  "nozzle  entrance"  is  not  critical  to  program  operation,  but  run  time 
is  reduced  if  a large  area  is  selected.  The  pressure  computed  by  proceeding  down 
the  port  must  agree  with  that  computed  (by  assuming  isentroplc  f;ow)  from  the 
throat  back  to  the  artificial  point,  within  a certain  tolerance. 

Values  of  c*  efficiency  and  throat  erosion,  which  are  output  from  the  bal- 
listics subroutines,  are  based  on  semi-empirical  correlations.  Thrust  may  be  cal- 
culated by  the  ballistics  subroutines  standing  alone,  but  requires  an  input  nozzle 
efficiency.  This  nozzle  efficiency  may  be  calculated  external  to  the  computer  por- 
gram,  using  such  empirical  expressions  as  contained  herein  or  elsewhere.  A 
series  of  motor  operating  parameters  and  running  integrals  are  output  and  saved 
for  use  in  other  modules.  If  the  other  modules  are  exercised,  the  thrust  parameters 
calculated  herein  are  overridden. 

Detailed  output  of  conditions  at  each  axial  node  in  the  calculations  may  be 
requested  as  a function  of  time  using  tabular  input  of  the  increments  desired. 

Calculation  of  the  bum-back  of  radial  slots  is  made  by: 

1)  entering  the  slot  faces  as  "nonbuming"  in  the  grain  design  program 
input; 

2)  overriding  the  grain  inner  diameter  calculation  with  the  outer  dia- 
meter in  the  slot;  and 

3)  explicitly  calculating  the  location  of  the  end  of  the  port  in  sub- 
routine MAIN3 

The  node  at  the  end  of  the  grain  retains  its  identity  number  but  is  then  a moving 
node.  As  the  burning  face  approaches  the  next  x-node  along  the  grain,  the  in- 
terior node  is  codec  out.  Table  input  is  provided  to  allow  for  ignition  delay  in  the 
slot  (radius  ignited  versus  time).  If  significant  delay  is  used,  the  face  of  the  grain 
will  not  remain  perpendicular  to  the  motor  axis  (axial  position  is  retained  at  the 
inner  diameter);  additional  table  input  is  then  needed  to  correct  the  grain  length 
as  the  port  bums  outward.  Evidence  that  corrections  are  needed  is  that  the  total 
propellant  weight  expended  will  be  under-predicted.  The  burning  rate  on  the  slot 
face  is  assumed  to  be  the  input  strand  rate. 

5.4.14  Subroutine  NER0DE 

This  routine  calculates  the  nozzle  throat  erosion  rate  from  input  tables 
or  from  a correlation  determined  from  analytical  calculations. 
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5,4,15  Subroutine  PRISM  1 


This  routine  calculates  the  initial  data  needed  for  each  specific  prism. 

Only  right  triangular  prisms  are  considered. 

5,4,16  Subroutine  PRISM2 

This  subroutine  has  been  written  for  use  with  the  Grain  Design  Program. 

The  purpose  of  the  routine  is  to  determine  the  two  points  of  intersection,  if  any, 
made  by  a vertical  line  and  a triangular  prism  located  in  three-dimensional  space. 

When  a prism  bums  out,  the  comers  bum  into  spheres,  the  edges  bum 
into  cylinders,  and  the  end  and  side  planes  bum  straight  out.  A line  will  inter- 
sect a prism  either  through  the  cylinder,  sections  (edges),  sphere  portions  (cor- 
ners), or  not  at  all. 

This  routine  considers  the  intersecting  line  as  always  a vertical  line  in 
the  standard  three-dimensional  system.  Hence,  the  vertical  line  is  defined  by 
an  X and  a Y vlaue.  The  number  of  line  intersections  necessary  depends  on  the 
particular  prism  beirg  considered  and  the  accuracy  desired  in  subsequent  computa- 
tions. 

Input  to  this  subroutine  consists  of  an  X value,  a Y value,  a T (bum)  dis- 
tance (positive  for  out-burning,  negative  for  in-burning),  and  those  computations 
made  in  PRISM1. 


5.4.17  Subroutine  PVS 

This  subroutine  is  used  only  with  the  ballistics  mode.  It  provides  peri- 
pheries, volumes,  flow  areas,  and  burning  surfaces  areas  (and  totals)  at  or  be- 
tween each  X-station  and  bum  interval  (time,  to  time  .). 

K rC  T 1 

Cross  sectional  flow  areas  at  timek  and  timek+1  are  supplied  to  this  ro- 
utine for  each  X-stat*on.  These  are  provided  in  two  tables  thus: 

Aj  j = Cross  sectional  flow  area  at  timek  and  X^ . 

A2  j = Cross  sectional  flow  area  at  timek+1  and  X^. 

The  n X-stations'  locations  are  provided  in  a table  of  X^. 
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Volume  burned  between  each  two  adjacent  X-stations  during  time,  to 

Jv 

timek+l  comPutec^  and  summed: 

Vtotal  = £ AVj 


AVj  f (A2,j  " Al,j*  + ^2,j+l  * Al,j+1^  ' ,Xj+l_Xj! 

The  average  periphery  of  the  flow  area  between  time^  and  time^+j , at 
Xj  is  calculated  as  the  change  in  flow  area  over  distance  burned  at  that  X-station. 
'Diis  is  calculated  only  to  be  printed. 


P = -U.  -1A 

J aBj 

The  average  burning  surface  area  during  the  interval  time^  to  timek+j, 
and  X^  to  X^+1  is  calculated  and  summed,,  For  this  calculation,  the  equation 
S=^V/^B  had  to  be  modified  because  is  not  necessarily  normal  to  S in  the 
ballistics  mode. 

5.4.18  Subroutine  READIN 

Subroutine  READIN  has  been  written  for  use  with  the  Grain  Design  Program. 
The  program  initializes  variables  and  reads  the  data.  A test  is  made  on  the  data 
in  order  to  pass  control  to  the  appropriate  place  for  storage  of  the  data. 

5.4.19  Subroutine  SFERE2 

This  subroutine  has  been  written  for  use  with  the  Grain  Design  Program. 
The  routine  will  determine  the  intersection  of  a line  with  a sphere  located  in 
three-dimensional  space. 

The  routine  solves  the  equation 

C = R2  - (Yc  - Y - (Xc  - X2)a 

where  R = The  radius  of  the  sphere. 

X^  = Denotes  the  location  of  a plane  in  three-dimensional 
space. 
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Y^,  = Denotes  a line  in  the  X plane. 

Y = The  Y coordinate  describing  the  center  of  the  sphere. 
Xj  = The  X coordinate  describing  the  center  of  the  sphere. 
If  C is  negative,  the  line  does  not  intersect  with  the  sphere. 

If  C is  zero,  the  line  is  tangent  to  the  sphere. 

If  c is  positive,  the  intersection  points  of  the  line  with  the  sphere  are  given  by: 

Z = z:  + C 

where  \ = The  Z coordinate  describing  the  center  of  the  sphere. 

5.4.20  Subroutine  STORE 

Subroutine  STORE  has  been  written  for  use  with  the  Grain  Design  Program. 
Its  purpose  is  to  determine  the  nature  of  the  input  data  and  store  it  in  locations 
to  which  the  main  program  has  access. 

Control  is  sent  to  the  program  from  Subroutine  READIN  by  the  statement 
CALL  STORE  (J,JO);  JO  has  been  assigned  an  integer  value  of  1 to  4 by  Subroutine 
READIN  depending  on  the  card  type  and  is  used  in  a "computed  go  to".  J is  a 
subscript  incremented  in  Subroutine  READIN  for  the  correct  storage  of  data  into  an 
array. 

5.4.21  Subroutine  TAB 

TAB  is  a routine  to  interpolate  (or  look  up)  values  from  the  tabular  inputs. 
Linear  interpolation  or  log-log  interpolation  is  used  and  no  extrapolation  is  made; 
rather,  the  value  at  the  boundary  of  the  table  is  returned.  The  tables  are  assumed 
of  the  form 

z = f (x,y) 

with  a call  statement  of  the  form 

CALL  TAB  (X,  Y,  Z,  No) 

where  No  is  a table  number  assigned  in  Subroutine  MAIN3.  The  size,  shape,  type 
of  interpolation,  and  order  of  table  loading  are  established  by  input  of  the  tables 
and  thus  are  not  transmitted  to  other  subroutines. 

Internally,  the  logic  is  established  to  search  the  table  by  stepping  from  the 
prior  point  to  the  new  one,  consistent  with  an  iterative  program  with  highly  re- 
petitious call  statements,  in  contrast  to  computing  indexes. 
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5,4.22  Subroutine  TABIN 

TABIN  is  a routine  which  prepares  input  of  tabular  data.  The  addressing  for 
the  tables  is  arranged,  indexes  are  initialized,  and  limits  to  table  size  established. 
Input  of  integers  and  non-tabular  program  inputs  are  simply  passed  back  to  READIN 
for  interpretation. 


5.4.23  Subroutine  UNI0N 

This  subroutine  has  been  written  for  specific  use  with  the  Grain  Design 
Program.  The  subroutine  performs  a union  with  Z points  solved  for  previously  and 
those  Z points  already  entered  in  the  table. 


A table  of  void  points  or  grain  points  is  kept  by  entering,  deleting,  or  in 
general,  performaing  a union  with  paris  of  Z points.  These  Z points  represent  the 
intersection  of  a line  and  the  basic  figures  allowed  for  in  the  Grain  Design  Pro- 
gram, Based  on  a flag,  the  Z points  can  be  treated  either  as  cuts  or  fills  to  the 
table.  Only  those  points  equal  to,  or  contained  within,  the  case  boundaries  are 
considered.  All  points  in  the  table  are  left  adjusted.  A count  of  the  number  of 
pairs  is  also  computed. 


Program  Symbols 
Fortran 


Description 


TABLE 


NOZE 

RDIN(2) 

RDIN(3) 

Z(l)  and  Z(2) 
ZEE  and  ZEET 
ZF 

J 

I 

I UFA 

KALI 

IDIF 

KUFA 


Location  where  void  or  grain  points  are  stored  after 
union  is  made. 

Number  or  paris  of  Z in  table. 

Lower  boundary  or  case  limit  points. 

Upper  boundary  or  case  limit  points. 

The  void  or  grain  points. 

Z(l)  and  Z (2) . 

Temporary  stroage  for  one  of  void  or  grain  points. 
Number  of  points  in  table. 

A count  of  number  of  points  already  looked  at  in  table. 

Flag  used  to  enter  or  delete  Z points.  This  flag  is 
also  set  outside  of  union  to  tell  if  table  contains 
grains  or  void  points. 

A counter  used  to  aid  in  entering  or  deleting  points 
from  the  table. 

A difference  value  used  to  left-adjust  points  in  the 
table. 

A branch  flag  used  within  the  program. 
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5.5  Link  40  Subroutines 


Overlay  40  contains  the  0DK  module  which  calculates  the  loss  due  to  finite 
rate  gas  phase  kinetics  as  incorporated  into  the  methodology  programmed  in  the 
SPP  code. 

This  overlay  includes  several  sub-overlays  which  are  executed  in  sequence, 
i.e.#  Link  42  is  executed  after  Link  41,  etc.  The  subroutine  descriptions  included 
in  this  section  are  in  alphabetical  order  within  the  individual  overlay  structure. 

All  of  the  subroutines  described  in  this  section  are  the  results  of  the  modi- 
fication of  the  0DK  program  described  in  Reference  3.  While  only  certain  of  the 
subprograms  of  that  reference  have  been  modified,  the  entire  pertinent  body  (i.e., 
0DK  documentation)  of  Reference  3 is  repeated  here. 


5.5.1  Subroutine  LINK40 

This  subroutine  consists  only  of  a call  to  subroutine  0DK  and  is  used  only 
to  facilitate  conversion  of  the  program  overlay  structure  to  other  computers. 


5.5.2  Subroutine  0DK 

This  subroutine  acts  as  the  driver  for  the  one  dimensional  kinetic  expansion 
calculation  (0DK).  It  calls  in  overAays  LINK41  through  LINK43.  Subroutine  0DK 
also  writes  (punches)  out  the  linkage  data  dealing  with  the  restricted  equilibrium 
option  of  the  0DE  module. 


5.5.3  Subroutine  STF 


This  subroutine  evaluates  the  thermodynamic  functions  Cp^/R,  H^/RT, 
S^/R,  fror  curve  fit  coefficients.  The  subroutine  uses  the  same  procedure  as 
subroutine  CPHS.  The  additional  functions  dfCp^/dT  and  free  energy,  G^/RT, 
are  also  computed.  The  calculated  functions  are  then  converted  to  the  internal 
computational  units  for  use  by  the  kinetic  expansion  calculations.  Also,  low 
temperature  thermod/namic  data  from  tables  is  used  when  required  if  the  LTCPHS 
directive  was  specified  in  the  input  to  the  SPP  code. 


5.5.4  Subroutine  STOICC 

For  each  reaction  this  subroutine  constructs  two  vectors  of  stoichiometric 
coefficients,  one  for  reactants  and  one  for  products.  Up  to  10  reactants  and  10 
products  may  be  considered  for  each  reaction.  The  total  number  of  entries  in  the 
resultant  linear  reaccion  table  is  600,  i.e.,  the  sum  of  all  stoichiometric  coef- 
ficients cannot  exceed  600. 


5.5.5  Program  Link  41 

This  subprogram  is  mainly  intended  to  facilitate  the  conversion  of  the 
SPP  code  overlay  structure  to  other  computer  systems.  Its  only  function  is  to 
call  subroutine  0DKINP  which  handles  the  input  to  the  0DK  module. 
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5,5/6  Subroutine  ECNV 


This  subroutine  translates  a BCD  string  of  characters  into  one 
floating  point  numeric  value.  E,I,  and  F formats  are  permitted  with  the 
result  always  a floating  point  number.  It  is  called  by  subroutine  REAXIN  to 
decode  numeric  fields  in  the  species  and  reactions  cards.  The  subroutine 
is  coded  entirely  in  FORTRAN.  A BCD  string  of  blanks  will  result  in  a floating 
point  zero  returned  value. 


5.5.7  Subroutine  NUMBR 

This  subroutine  converts  a one  character  BCD  number  to  a FORTRAN  integer 
number.  The  subroutine  is  coded  entirely  in  F0RTRAN. 
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5.5.8  Subroutine  CfoKINP 


This  subroutine  provides  the  input  processing  for  the  kinetic  expan- 
sion calculation,  ft  performs  the  following  functions: 


1)  Variable  Initialization  to  nominal  values 

2)  Calls  subroutine  REAXIN  to  input  the  reactions  cards  and  species 
cards  if  necessary 

3)  For  an  0DE-0DK  problem,  calls  subroutine  SELECT  to  select 
those  species  to  be  considered  for  the  kinetic  expansion 
calculation 

4)  Reads  $0DK  namelist  input  data 

5)  Converts  nozzle  geometric  parameters  from  input  units:  inches, 
degrees;  to  internal  computational  units:  feet,  radians 

6)  Computes  nozzle  tangent  coordinates  using: 


rt  = 1 + Rd  (1  - cos  6) 
xt  = Rd  sin  0 


7)  For  conical  nozzles,  computes  the  axial  coordinate  for  the  exit 
station  from  the  following  relation: 


xexit 

xexit 


rt  + xt  • tan  6 
tan  0 


- ( 1 + Rd  -VEP) 


2 


xexit 

xexit 


<XA 


8)  ' For  conical  nozzles,  the  internal  axial  print  stations 
are  computed  using: 


X,  = 


VaRPRNT(J)  - rt  + xt  • tan  6 
tan  0 


yj*xt 


X,  = 


R 2 - ’ 

d 


1 + Rd  - (ARPRNT(J) ) 


4 


Xj  <Xt 
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9)  The  sum  of  input  or  selected  species  concentrations  is  checked 
for  unity  ( + XMFTST,  where  XMFTST  is  an  input  number),  and 
then  normalized. 

10)  11  the  input  parameter  RZN0RM  is  input,  the  input  contoured 
nozzle  table  is  normalized  by  RZN0RM. 

In  addition  to  the  above,  the  following  modifications  have  been  made  to 
the  0DKINP  routine. 

a)  If  the  nozzle  geometry  has  been  input  via  the  $GE0M  namelist,  0DKENP 
traps  and  stores  the  variables  from  this  input  into  the  correct  locations 
and  performs  the  appropriate  unit  conversions. 

b)  Was  modified  to  accept  the  parabolic  and  circular  arc  options  (IWALL=2 
and  3)  from  the  $GE0M  namelist.  This  was  done  by  computing  the 
nozzle  wall  points  and  derivatives  at  20  equally  spaced  points  and 
substituting  them  into  the  tables  used  in  the  contoured  wall  input 
option. 

c)  Was  modified  to  accept  as  input  the  quantities  needed  to  calculate 
the  zero  particle  lag  flow  equations. 
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5,5,9  Subroutine  REAXIN 


This  subroutine  processes  SPECIES,  REACTIONS,  and  THIRD  B0DY  REAX 
RATE  RATT0S  input  cards.  Reference  may  be  made  to  Volume  III,  Section  2.3.2, 
the  Program  User's  Manual,  for  a complete  description  of  input  requirements.  A 
table  of  all  species  appearing  in  the  input  reaction  set  is  generated  for  further 
processing  by  subroutine  SELECT  if  required. 


5,5,10  Subroutine  SELECT 


This  subroutine  provides  the  interface  logic  required  to  select  the 
minimum  species  list  required  for  the  kinetic  expansion  calculations.  The 
subroutine  is  only  used  for  the  0DE-0DX.  interface.  The  list  of  all  species 
appearing  in  the  input  reaction  set  is  matched  against  the  list  of  species  con- 
sidered for  the  equilibrium  calculation.  All  species  which  appear  in  both  a 
reaction  and  the  equilibrium  calculation  list  are  selected  for  the  kinetic  expan- 
sion calculation.  If  a species  appears  in  the  reaction  set  but  has  not  been 
considered  for  the  equilibrium  calculation,  the  program  prints  an  error  message 
and  terminates  the  current  case. 

If  the  INERTS  directive  was  specified  those  species  specified  under 
that  directive  will  be  added  to  the  list  for  the  kinetic  expansion  calculation  along 
with  the  other  species  selected. 

If  the  INERTS  directive  was  not  specified,  all  those  species,  con- 
sidered for  the  equilibrium  calculation,  whose  mole  fractions  are  greater  than 
or  equal  to  an  input  selection  criterion  will  also  be  selected  for  the  kinetic 
expansion  calculation.  Species  selected  in  this  way  will  be  listed  as  inert 
species  on  the  program  output  since  they  do  not  enter  into  chemical  reaction. 

This  routine  was  modified  to  select  pairs  of  condensed  phases  if 
either  of  the  phases  passes  one  of  the  above  selection  criteria. 
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This  subprogram  is  mainly  intended  to  facilitate  the  conversion  of  the  SPP 
code  overlay  structure  to  other  computer  systems.  Its  only  function  is  to  call 
subroutine  PACK.  See  Section  5.5.13  below. 
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5.5.12  Subroutine  CQfNVRT 


This  subroutine  converts  input  data  from  the  externally  input  units 
to  internally  used  computation  units.  In  order  to  conserve  computation  time 
during  the  kinetic  expansion,  parameters  such  as  molecular  weights,  are  in- 
cluded in  these  conversions.  Primed  numbers  are  input  quantities. 


a)  Reaction  rate  ratio  input  for  reactions  requiring  third  body  terms 
units:  unities s 

internal  units:  ( lb  s-mas  s/lb -mole)"1 
formula:  XMM.  . = XMM.  ./Mw 

*■ 

b)  Pre-exponential  reaction  rate  parameter 
input  units:  cm,  °K,  g-mole,  sec 

Q 

internal  units:  ft  , 0 R,  lb-mole,  sec 


aj = 


aj  * (.0160183)**  1.8nj 
n 

n 

i=i 


n MV/ij 


Where  X depends  on  the  order  of  the  reaction. 


and  .0160183  = 


3.531  • 10'5  ft3 

, 3 

1 cm 


1 g-mass 


2.2*  10  3 lbs-mass 


c)  Exponential  Term: 

input  units:  kcaJ/mole 
internal  units:  °R 
formula:  b.  =bj  • 905.770 


«*iiere  905.770 


1000  cal 
1 kcal 


1 I 8°R 

1.98726  cal/mole-:K  TTIFk 
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<9  Equilibrium  Constant  Multiplicative  Factor: 
Input  units:  not  input 

o 

internal  units:  (Ibs-mass)  - °R/ft 
formula: 

n 

i-i Mw. Vii 

DkTEF(J)  = —i-i 

n 

II  0.73034 


where  .73034  = 49,721.011-  A -goun.d\ ^ 

(lbs-mole)  R 


1 atmos  2 

68,059.59  poundals/ft^ 


e)  Pressure: 

input  units:  PSIA 

9 

internal  units;  poundals/ft 
formula:  P = P’  - 4633.056 
where 

2 

4633.056  = -14-4  *n  • 32.  174  -&—0 
1 ftz  sec 


f) 


The  initia1  reference  enthalpy  is  computed  using 


i=l 


Subroutine  C0NVRT  has  been  modified  to  calculate  the  total  amount  of  con  • 
densed  phase  present  at  the  kinetic  expansion  initial  conditions. 
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$>5,13  Subroutine  PACK 


On  the  basis  of  those  species  currently  being  considered,  this  sub- 
routine packs  species  and  reaction  information  from  the  master  tables  into 
those  control  sections  utilized  by  the  one-dimensional  kinetic  expansion  sub- 
program. 

The  following  is  a sequential  description  of  the  packing  procedures: 

1)  Thermodynamic  data  for  the  species  being  considered  is 
reed  into  core  storage. 

2)  The  chemical  species'molecular  weights  are  computed 

3)  The  symbolic  reactions  are  checked  for  mass  balance. 

4)  For  a contoured  nozzle  the  slope  at  each  input  wall  point  is 
computed  using  subroutine  SLP.  The  wall  coordinates,  and 
each  computed  slope  are  printed  for  each  input  wall  point  and 
the  print  stations  are  set  to  the  input  axial  coordinates. 

In  addition  to  the  above,  the  PACK  routine  processes  condensed 
phases  by  setting  R^O  for  those  species  and  storing  pointers  up  to  10  condensed 
phase  species  (5  pairs).  This  routine  has  also  been  modified  to  calculate  the  print 
positions  for  specified  area  ratios  when  the  contoured  (spline)  wall  option  is 
selected. 
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5.5.14  Subroutine  PRES 


This  subroutine  is  used  (when  JPFLAG  s 1)  to  compute  the  derivatives 
of  an  input  pressure  table. 

This  subroutine  is  also  used  (when  JPFLAG  = 0)  to  generate  a pressure 
table  through  use  of  an  average  expansion  coefficient,  Ne.  The  generated  table 
extends  from  the  initial  contraction  ratio  through  the  nozzle  attachment  point  plus 
one  normalized  throat  radius. 

Input  Pressure  Table  Derivative  Computation  (TPFLAG  = 1) 

If  a pressure  table  of  NTB  entries  is  input,  the  table  of  first  deriva- 
tives is  computed  using: 


, 1 < n < NTB 
, n « NTB 

The  pressure  at  the  initial  axial  position  is  obtained  by  interpolation 
using  subroutine  SPLN. 

i 

i 


dP 

dx 


dP 

dx 


x 


h 


XNTB 


: Vi> 

x(n+l)~  X(n-1) 
x(n)  " X(r«-1) 
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Internally  Computed  Pressure  Table  Computation  (IPFLAG  = 0) 


An  average  equilibrium  pressure  expansion  coefficient  from  the 
chamber  to  the  throat,  is  commited  by  iteration  using  subroutine 
SUBNE.  The  initial  value  for  N (1'  is  1.2. 


The  approximate  equilibrium  contraction  ratio  at  the  initial  axial 
position  is  computed  from: 


where 


= pressure  at  the  initial  axial  position 

P = equilibrium  chamber  pressure 
c 


A check  is  then  made  to  determine  the  compatibility  between  the 
nozzle  geometry  and  the  requested  contraction  ratio. 


ifcT* 1 +[Ru+RiHi  -cos  %]• 

the  circular  arcs  Ru  and  R^  overlap  and  the  following  error  message  is  printed: 

INLET  GEOMETRY  INCOMPATIBLE  WITH  INITIAL  CONDITIONS 
The  program  will  proceed  to  the  next  case. 
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Tables  for  pressure  and  its  derivatives  are  constructed  as  functions  of  area 

ratio,  a,  and  expansion  coefficient,  Ne.  Formula 

used  for  the  j + l iteration  for  pressure  is:  • * 


V1 


P0+1)  p(l> 


fi-rfi  ■ 


-1 


N -1 
e 


[2/(N  +1)] 


(Ne+l)/CNe-l) 


2/N 
(J)'  6 


1 - 


<Ne-l)/Ne 


-1/2 


a -i 


for  J = 1 


,(1): 


n+l 


p_!  d(P/Pc  ) 

Pc  I + 

'x 

n 


(x 


n+l 


V 


where  n refers  to  the  n ^ table  entry. 


The  pressure  derivative  formula  used  is: 


d(F/Pc) 

dx 
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.7/>~ 


.r 


Next,  tables  for  pressure  and  its  derivatives  are  constructed 
by  the  program.  Table  entries  are  at  increments  of 


-x/7S 

for 

Xj  < x < 0 

(Rd  sin  e)/25 

for 

0 <x  < R^  sin  0 

and 


1/25 


for 


R^  sin  0 < x < sin  p + 1 


where  the  initial  nozzle  axial  position,  x1#  is  computed  from: 


- 

•\fei~-  1 - (R  + R=)  • (1  - cos  e.) 

ci  = " (Ru  + Ri)  ‘ sln  6i  + teTeJ 


See  Figure  5.5-1  below: 


Area  ratio  and  its  derivative  and  (a  and  are  found  by  the  five  formulae 

below: 


1)  x < Xj  + Rj  sin  0^ 

js. . ~ 2 lx~** 1 .vr 

«* 


2)  Xj  + Rj  sin  0j  < x < - Ru  sin  Qj 


a = “ Ri^  _ cos  “ ^x_xi  “ Ri  sin  ei^  tan  6i] 

~3x”  = “ 2 *>/«"•  tan0i 


3)  -Ru  • sin  0 i < x < 0 


- [■  * r»  (>  -f-  -%) 


4)  0 < x < * sin  0 


da  _ 2x  <JT 
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5)  Rd  • sin  0 < x ^ Rd  • sin  e + 1 
for  cone 

a * £rt  + (x-x^  * tan  0 

jSS-s  2 • a • tan  0 
dx 


for  contour 


a ®Y 


2 


da  _ 
dx 


2*  Y* 


dY 

dx 


Three  special  points  are  included  in  the  pressure  table.  These  are 
a point  at  x = x4  such  that 


d(P/Pc) 

dx 


and  two  points  at  x = 0 such  that 


d(P/Pc) 

dx 


equilibrium 

N 

e_  . 

•ys* 


3Ne  ~ 1 

2(N^ry 

r,  2,  i 

LNe+1  J 


with  R*  = Ru  and  R*  = Rd,  respectively. 
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The  following  Items  are  Input  directly  to  the  computer  program  as  described 
in  Volume  III,  Section  2.8.5,  Figure  2-2,  left  to  right. 


u 


RI 

RWTU 

THETAI 

THETA 

RWTD 
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The  purpose  of  this  subroutine  is  to  supply  derivatives  for  a tabulated 
function.  The  end  point  derivatives  may  be  specified  or  are  calculated  inter- 
nally by  parabolic  interpolation.  Interior  point  derivatives  may  be  found  by 
a cubic  spline  fit  procedure. 

Calling  Sequence: 

X is  a table  of  Independent  variables,  x 

Y is  a table  of  the  dependent  variables , y^ 

N is  the  number  of  entries  in  each  of  the  tables  X,  Y,  and  YP. 
1=  1,  ...N 

MFLAG  this  entry  is  a flag,  m,  such  that 

m > 0 Implies  x is  equally  spaced 

m < 0 implies  x is  not  equally  spaced 

|m|  = 1 y'  will  be  continuous 
|m|  = 2 y*  arid  yM  will  be  continuous 
YP  is  a table  of  the  derivative,  yj 
W1  working  storage  of  length  N 
W2  working  storage  of  length  N 


below: 


W3  working  storage  of  length  N 


IFLAG 


this  entry  is  a flag,  i,  such  that 

i = o implies  value  for  YP(1)  and  YP(N)  will  be  calculated 
internally  by  parabolic  differencing 

i = 1 implies  values  for  YP(1)  and  YP(N)  will  be  input 


i=  1 


Method 


The  cubic  spline  fit  procedure  utilizes  the  interpolation  formula  given 


3 2 

y - A(x  - x ) + B(x  - x ) + C(x  - x ) + D 
7 ' o £ o o 

y*  = 3A(x  - xq)  + 2B(x  - xq)  + C 
yn  = 6A(x  - xq)  + 2B 

The  piecewise  cubic  fit  to  a tabular  function  by  the  above  relations  will  yield 
a discontinuity  in  the  second  derivative  yM,  between  adjacent  fits  of: 


1 


tion  equal  to  zero  so  that  the  second  derivative  is  continuous  across  Juncture 
points.  As  applied  to  a tabular  function,  the  above  procedure  results  in  a set 
of  linear  simultaneous  equations  (tri-diagonal)  to  be  solved  for  the  y| , provided 
that  values  for  y'  at  the  end  points  are  known. 
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5,5,16  Subroutine  SUBNE 

Calculates  the  average  equilibrium  pressure  expansion  coefficient 

from  the  chamber  to  the  throat  by  iteration  from  the  following  formula 
(Newton's  method): 


where  -1.2. 

P*  is  the  equilibrium  throat  pressure 
Pc  is  the  equilibrium  chamber  pressure 

This. subroutine  Is  used  by  subroutine  PRES. 
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5,5.17  Program  Link  43 


This  subprogram  is  mainly  intended  to  facilitate  the  conversion  of  the  SPP 
code  overlay  structure  to  other  computer  systems.  Its  only  function  is  to  call 
subroutine  MAIN  ID  which  handles  the  control  of  the  integration  of  equations  de- 
scribing one  dimensional  reacting  gas  flow  through  a rocket  nozzle. 


5.5.18  Subroutine  EF 

This  subroutine  computes  equilibrium,  constants,  K 


j 


« 


Kj  EK(J) 
also  computed  are 


DATEP(I) 
'p  A 


dK. 

l 

dT 

DKT(J) 

where: 

Fti 

Hti 

DATEF  (J) 


l‘U  + & R, 


n 

E 


Ht 


i 


= species  free  energy  at  the  current  temperature 
= species  enthalpy  at  the  current  temperature 
= is  defined  in  subroutine  C0NVRT,  Section  5.5.12. 
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5.5,19  Subroutine  DERIV 

This  subroutine  computes  the  total  derivatives  and  the  partial 
derivatives  g.  ^ described  in  the  analysis  presented  in  Section  3 of  Reference 
3* 

The  implicit  integration  method  used  to  integrate  the  differential 
equations  governing  the  chemical  system,  i.e. 

yi  = Vx,yl*  ' ' *yNSP+3*  i=1'  ' * " NSP+3 

v/here  the  variables . 


yt  i=l,  . . NSP+3 

are  V,  p,  T,  i=l,  . . NSP 

respectively,  requires  evaluation  of  the  Jacobian  of  the  system,  i.  e. 


Bf 


pij”  ay 


i 


i=l,  . . .,  NSP+3 
j=l,  . . NSP+3 


Subroutine  DERIV  computes  only  certain  of  the  p (those  taken 
with  respect  to  C^)  and  the  others  are  computed  in  subroutine  FLU. 

Also  calculated  by  DERIV  are  the  reaction  rates,  kj,  and  the  net 
production  rates,  X^. 

The  generalized  chemical  reaction  which  is  handled  by  this  sub- 
routine is  defined  by: 

NSP  NSP 

. ,?!  nj  V 

where  Mj  represents  the  ith  chemical  species. 

The  reverse  reaction  rate  constant  is  defined  by  the  equation: 


kj  SK(J) 


=aJ  • T ~n.j  ■ exp  (-bj/T) 


The  net  production  rate  for  a reaction  is  given  by: 


Xj  X(J) 


f NSP  u NSP  „•  I 

iv  w -- -H,ci  'J 
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where:  X depends  on  the  order  of  the  reaction* 


NSP 

with 

MJ  = 1=1  mMi.i  ' C1 

for  reactions  requiring  a third  body 

and 

Mj  = 1 

for  all  other  reactions 

The  net  individual  species  production  rate  is  given  by  the  equation: 
dC.  _ ^ 

-sr  wo  -v  fa 


where: 

Ka  = (Mw1  • p • r*)/V 

fyEVvij 

The  partial  derivatives  of  the  net  species  production  rate  with 
respect  to:  the  chemical  species;  the  gas  velocity;  the  gas  density;  and  the  * 
gas  temperature  are:  . 

_ T »x, 

BT(I,K)  = Kj  • i = l NSP 

1 k = 1,  . . NSP 

1 dCi 

PHI(1, 1)  = -“7"-^-  i = l NSP 

i dC.  _ y 3X 

PHKI.2)  = ■ + j=i  ^1  = 1 NSP 

l 

_ y ax 

PHI(1, 3)  = Ki  j=l  ST  i NSP 


* X = (Vjj  - v^)  so  that  X = 0 for  binary  exchange,  X = 1 for  most 


i=l 

dissociation  recombination  reactions. 
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The  subscript  notation  used  above  is: 

1 *=  Species  subscript 

j « Reaction  subscript 

i = Total  number  of  chemical  reactions 

m * Number  of  reactions  requiring  third  body  terms 

NSP  = Total  number  of  gaseous  species 

Subroutine  DERIV  calls  subroutine  FLU  to*  calculate  the  derivatives  and 
partial  derivatives  of  V,  P,  and  T.  In  the  event  of  solidification , this  routine 
also  recalculates  the  condensed  phase  species  derivatives  and  their  partial  de- 
rivatives and  recalls  subroutine  FLU  to  recalculate  the  flow  derivatives.  See 
Volume  I,  Section  for  the  description  of  this  procedure  and  the  equations 

calculated. 
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5,5.20  Subroutine  FLU 


This  subroutine  computes  the  total  derivatives  and  the  partial 
derivatives  a 4 and  3^  for  the  fluid  dynamic  equations.  While  the  flow  is 
subsonic,  pressure  defined  fluid  dynamic  equations  are  used.  When  the 
flow  becomes  supersonic,  area  defined  fluid  dynamic  equations  are  used. 

The  summation  terms,  energy  exchange  term  B,  the  diabatic  heat  addition  term 
A,  the  Mach  number,  and  all  the  partial  derivatives  of  these  terms  are  com- 
puted. For  a subsonic  integration  the  pressure  and  its  derivatives  are  obtained 
from  the  subsonic  pressure  table.  For  a supersonic  integration, the  area  ratio 
and  its  derivatives  are  computed  from  the  input  geometric  constraints. 

The  calculations  logically  fall  into  three  types:  a)  Those  done  for 
all  integrations;  b)  Those  done  only  for  subsonic  integration;  c)  Those  done 
only  for  supersonic  integration.  The  following  will  adhere  as  closely  as 
possible  to  a sequential  description  of  the  computations. 

The  operators  $ (i , 1)  , 1 = 1,2,3  are  defined  as 
$ (i 1 1)  « P (C.,  V) 

ft(i,2)  = )3  (CA  , p) 

$(i.3)  = fi  (Cj.T) 

dC{ 

The  total  derivatives,  f i = -gj-  , for  i = 1,  . . . , n are  computed  as 
tt>4r* 


where  n is  the  total  number  of  species,  NSP. 


Computation  of  the  Summation  Terms  and  their  derivatives; 
First  Summation 


81 

SI 

- 

1 11  C 
I -S  - 
R 1=1 

dSl 

dV 

DS1V 

a 

1 11 

B.  * 

dSl 

dp 

DS1R0 

St 

1 11 
rSi 

dSl 

dT 

DS1T 

a 

1 11 

tF-i'i 

dSl 

dCi 

- DS1C{I) 

- 

i . fs 

R Lj=l 

Second  Summation 

, n 

82 

02 

B 

-i-s 
r-t  i=: 

dS2 

dV 

DS2V 

S 

i n 

-i-  -S 
R-T  1=1 

dS2 

dp 

DS2R0 

= 

i n 

r • t i=; 

dS2 

dT 

DS2T 

a 

_j_.£ 

R-TW 

'1 


R. 


ffr  DS2C(I) 


R •[|5l,’(C).01)  VS1'Kl]  " 


dc. 


1 £ *-  , — 1^  1 S2 

R-TiM  t_*(1.3) ' hl  + dx  Cpl  j T 


1 


^ n.  ^(C  ,C  )'hl 

R-Lf-l  I » 


*R  is  the  gas  constant/molocular  wt.  of  species  i except  for  condensed  phases 
whfcre  R^O. 
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Computation  of  the  Energy  Exchange  Term  B and  its  Derivatives: 


B 

BB 

- 

1=1  . 

7 

S2 

dB 

dv 

DBBY 

- 

1=1  . 

7 

as2 

av 

dB 

dp 

DBBR0 

MS 

1=1  . 

7 

as2, 

ap 

dB 

dT 

DBBT 

= 

1=1  . 

7 

as2 

ai 

. J2  . & 

y2  aT 

dB 

dCi 

DBBC(I) 

= 

1=1  . 
7 

as2 

act 

+ S2  . a* 
y2  dci 

Computation  of  the  Diabatic 

Heat  Addition  Term  A 

and  its  Derivatives 

A 

AA 

SET 

Sl-B 

dA 

DAAV 

asi 

dB 

dV 

SS 

av 

dV 

• 

dA 

dp 

DAAR0 

= 

asi  _ 
8p 

dB. 

dp 

dA 

dT 

DAAT 

asi 

dB 

SS 

St 

dT 

‘ 

dA 

d^ 

DAAC(I) 

= 

asi 

3 ci  " 

dB 

dci 

i= 

I 
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' 3 


m 


Computation  of  the  Mach  number  and  its  derivatives: 


M2 

XM2 

II 

»-3 

*M2 

av 

DM2V 

_ 2 -M2 

V 

3m2 

3T 

DM2T 

2 2 
Ml  _ MT 
T y 

±7 

3T 

3M2 

*C| 

DM2C© 

XKi  r ay  i 
- - M Lac;  r 

Rin 

+ TJ 

i **  1, . , 

> • , n 

For  the  subsonic  portion  of  the  nozzle,  pressure  defined  fluid  dynamic 
equations  are  used.  The  pressure,  and  its  first  and  second  derivatives  are 
computed  via  interpolation  in  the  pressure  table  generated  by  subroutine  PRES. 
The  Subsonic  Gas  Velocity  derivatives  are  computed: 


dV 

dx 

FNX(l) 

1 

sr  — ■ • 

P‘V 

dP 

dx 

d[FNX(l)] 

ax 

AL(1) 

1 

“ " 'p  V * 

d2P 

dx2 

0(V.V) 

BETA  (1,1) 

1_ 

V 

dV 

dx 

U(V,p) 

BETA(1,2) 

i 

p 

dV 

dx 

The  Subsonic  Gas  Density  derivatives  are  computed: 


dx 

afrNXfe)] 

dx 


FNX(2) 

AL(2) 
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fi{ P.V)  BETA(2 , 1)  p • av 

/J{p,p)  BETA(2,2)  ■-  ^ ~ p ' bf 

P(p,T)  BETA (2, 3)  =-  P • f*  ~ p^2‘ fj  * £ 

% 

Wp.Cj)  BETA (2 , 1+3)  = - “2^  ' ^ * JJ-  p-|*'  1=1..  ...n 
The  Subsonic  Gas  Temperature  derivatives  are  computed: 


dT 

dx 

FNX(3) 

T • 

SLngBU 

dx 

AL(3) 

XzL  . 
y 

/KT.V) 

BETA(3 ,1)  = - 

T • 

0(T.p) 

BETA(3 ,2)  = - 

T • 

/J(T.T) 

BETA(3,3)  = 

I . IT 
T ’ dx 

^(T.Cp 

BETA(3 , 1+3)= 

T • . 

~Z±  . i . dP  _ -> 
L y P dx  J 


1 !" d2P  /flP  N2  I 1 

p ' Ldx2  ’ P J 

aj 

av 


dp 


+ T 


_i_  d£ 
y2.p  dx  8T 


51 

BT 


r JL  dp  5z  M 
Ly2.p  dx  dCA- 
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For  the  supersonic  portion  of  the  nozzle,  area  defined  fluid  dynamic 
equations  are  used.  The  area  ratio,  and  Its  derivatives  are  computed  according 
to  the  Input  geometric  constraints. 

Area  ratio  and  Its  derivatives: 

1)  On  the  circular  arc  of  radius  (input  Item  RWTD)  defining 
the  downstream  throat  region,  X s ^tangent 


2)  For  a conical  nozzle  and  X>  X 


tangent 


da 

dx 


j't  + (x  " xc)  tan  °t  J 

. 2 + ^x-x^  tan  °t  j'tan  ®t 


4"  2 tan^  0 

dx*  S 
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3)  For  contoured,  parabolic,  and  circular  arc  nozzle  and  x>^angent 


a - Y2 


where  Y,  dY/dx  are  computed  via  interpolation  in  the  tables  of  Y and  its  deriva- 
tives generated  by  exact  calculations  or  by  subroutine  SLP. 

The  Supersonic  Gas  Velocity  derivatives  are  computed: 


dV 

dx 

FNX(l) 

JL  . ['I  da  _A“1 
m2-i  La  J 

ilm xjiJl 

dx 

AL(1) 

V_  . i . rA  _l.p  da. 
M2-l  a L dx2  8 V 48 

■)’] 

/»(V,V) 

BETA(l.l)  = 

I . dV  JL  . dV  . dM? 

V dx  ^ ^ dx  dV 

V . dA 

■ 

2 by 

M -1 

W.  p) 

BETA(1 , 2)  = 

V &A 

"m2-i’  ap 

BETA(1 ,3)  = 

1 dV  . dM2  V . 

" m2-i  ■ **  aT  m2-i 

dA_ 

dT 

BETA(1 ,1+3)  = 

1 . dV  . dM2  JL  . 

“ m2-i  aci  ’ m2-i 

dA 

dCj  W" 
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Ihe  Supersonic  Gas  Density  derivatives  are  computed: 

. FKX(2)  » - p{  ”A)  + A] 


ttFjTOJ  ^(2) 


M 

P • -y~ 
M -1 


1 f da  I /"da^  l 

a L j 2 a Vlx^  j 


0(p,V) 

BETA(2,1)  = 

r 1 f 1 da  A 3M.  _J_  . \h.  "j 

> V.,l>  ■ = -V  *>  Mz-1  wJ 

/J(  P/p) 

BETA(2  / 2)  = 

JL.de.  + _fi_  . |A 
p Mz-1  p 

Mp.T)  ' 

BETA(2,3)  = 

r 1 ,M  da  y&M2  1 • — 1 

P<L(m2-i)2  dx  ^ dT  m2-i  9tJ 

PtP,^) 

BSTA(2,l+3)  = 

- 1 /■  1 da  -\dM2  , _i .^A  ~j 

11  ' i- (M2-!)2  ^ * dX  790‘  mW 

l r l / • • • / 

The  Supersonic  Gas  Temperature  derivatives  are  computed: 

FNX{3)  - - T ' [(r-D  ' -^T  v a ” kJ  + B j 
dx  L M -1 


mmi  alO, 


S(T,V)  BETA(3 , 1)  = T 


0(T,p)  BETA(3 , 2) 


2 , _ ,2 
M y-1  ,!  d_a. 

M2-l  8 dx2 

_I.f  d2.Y  I 
a V dx  ^ J 

r y - 1 V 1 da  A 

»v  *r  V-l 

. §A 

a_B*i 

V-D2Vadx  ' 

dV 

“av  J 

, r i m2  .M 

aj  i 

v1  2,  dp 
L M -1 

ap  j 
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P(T,T)  BETA(3,3) 


I.M  7.1*4  . — lli. 


r y-i  r Ida  -s  sm2 
L^zUdx  " J 3T 


+r 


, , MZ  3A  93 
-1  • — r*  • »r  - sr 


M2rl 


9T  91 


/3(T,  Cj)  BETA  P , 1+3)= 


_ mL  . ri  da 

„.2  \a  dx 
(M  -1) 


T.r  * 

L (M2-!)2  ^ “ * 


M2  , 3A 


9B  f\  da  a>iZ_1 


1 = 1, 


• • • / D 
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v: 


5 >5 .21  Subroutine  GTF 

This  subroutine  computes  the  effective  gas  constant,  gaseous  heat 
capacity,  y,  dy/dT,  dy/&  from  the  following  formulae: 

NSP 

R - 2 C.  • R 

1-1  1 1 

NSP 

Cp  - S c.cp. 

1 = 1 1 1 . 

v , £J2_ 

7 Cp-R 


221 


« - 2L.&1-U 

Cp 


NSP 
s CL 
1=1 


acpj 

Tt" 


221 

sc 


y.  fr-1) 


Cp  J 


i - 1 , • • • , n 


For  condensed  phases,  R^O,  to  account  for  the  assumption  that  the  particles 
exert  no  pressure  on  the  gas. 
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5.5.22  Subroutine  IAUX  (1IL,  H,  QK,  RK,  IX) 

This  subroutine  performs  implicit  integration  according  to  the  method 
discussed  in  Section  of  Reference  3.  The  Increments  for  the  chemical  species 
concentrations  and  the  fluid  dynamic  variables  at  thelorward  point  are  cal- 
culated by  solving  the  appropriate  implicit  finite  difference  formulas.  Subroutine 
IAUX  also  performs  explicit  integration,  using  a modified  Euler  method,  when 
the  gas  temperature  falls  below  an  input  value. 

The  calling  sequence  parameters  are: 

HL  - last  integration  step  size 

H - current  integration  step  size 

QK  - last  increments  for  variables 

RK  - computed  increments  for  variables 

JK  - 1 initial  3 steps 

2 general  step 

3 special  step 

4 restart  step 

The  total  derivatives,  f.  , and  partial  derivatives,  p.  . at  the 

i,n  i , J , n 

back  point  are  calculated  in  subroutines  DERIV  and  FLU. 

The  special  step  calculation  is  used  at  print  stations,  in  halving 
the  step  size  if  required,  or  for  integrating  to  specific  calculation  stations.  If 
the  special  step  calculation  is  used  to  determine  the  properties  at  a print  station, 
the  calculation  is  resumed  using  the  general  step  calculation  and  the  previous 
step  size. 

After  each  integration  step,  subroutine  IAUX  obtains  the  derivatives 
at  the  then  current  axial  position. 


5-94 


For  Implicit  Integration  the  equations  used  are: 
Initial  Step  and  Restart 


°i,oh 


N 

+ r 


e.  . K 

i.j.o  j 


General  Step 

^ifn+l 

Special  Step 


&..* 2 1 


f . + a.  h 

i,n  i,n 


N 

+ r 


i,j,n  J ,n+l 


n+1 


i,n+i  (2h  +h  )-h 
ui'JL  n n 


N 


k . + 

i,n 


f . + a.  h 

i,n  itn  n+1 


» T n 

4-  r A 1 r - n / 

J-l  1>j>n  j’n+1]  hn+l  ' 


h..+h  ) 
n+1  n 


For  explicit  integration  the  above  equations  are  used  deleting  the 
partial  derivative  terms  a and  £ . 

If  the  option  to  generate  input  tables  for  the  Turbulent  Boundary  Layer 

Nozzle  Analysis  Computer  Program  was  selected,  tables  of  M,  P/P  T/T  , C , 

c c p 

V,  p,  are  tabulated  using  subroutine  TABGEN.  While  the  code  remains  in  this 
routine  to  do  the  above,  the  SPP  code  does  not  use  these  values  in  any  way. 
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5.5,23  Subroutine  INT 


Provides  control  for  the  implicit  Integration  procedure , determines 
the  proper  set  of  nonhomogeneous  equations  to  solve,  and,  after  each 
integration  step,  computes  the  next  integration  step  size  according  to 
the  following  relations: 


w - 2W 


hn+2  * 2 n+1* 


kt,n41  ~ 2kj.n  4 kl.n-l 
3k,  ....  - k. 


i.nfl  i,n 


< 

I MAX 


A»n+1  i.n  l.n-1 


3k  - k 
1,1*1  Ki,n 


■MAX 


k - 2k,  + k.  , 

_ i,n  i.n-1 


3k  - k, 
i,n+l  i,n 


'MAX 


6_ 

10 


6 


6 


On  option,  (JF=1)  only  the  fluid  dynamic  variables  are  used  in 
determining  the  next  integration  step  size. 

If  the  step  size  is  halved  for  the  fourth  step,  the  integration  is  re- 
started using  one-half  the  original  step  size. 

The  correspondence  between  equation  number  and  physical  property 
is: 

Equation  Number  Property 

1 Velocity  of  Gas 

2 Density  of  Gas 

3 Temperature  of  Gas 

4 *♦  NSP+3  Gaseous  species  mass  fraction 

(1  -*NSP)  corresponds  to  (4  -♦  NSP  + 3) 


When  the  flow  is  supersonic,  continuity  is  used  to  control  the  integration  step 
size  to  insure  that: 


C pVA)  " ( pVA)^ 
<-PVA)n+1 


< CtfNDEL 


where  CGfrJDEL  Is  an  Input  relative  criterion. 
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5, 5 >24  Subroutine  LESK  (Y) 

This  subroutine  Is  a single  precision  linear  equation  solver  which  is 
used  to  perform  the  matrix  inversions  required  by  subroutine  IAUX.  Gaussian 
elimination  is  used  with  row  interchange  taking  place  to  position  maximum 
pivot  elements  after  the  rows  are  initially  scaled. 


l 

s 


l 

T 
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5,5.25  Subroutine  MAINin 


This  subroutine  provides  the  overall  logic  control  for  the  one-dimen- 
sional kinetic  expansion.  The  following  functions  are  controlled: 

1)  Variable  initialization 

2)  Option  to  start  the  kinetic  expansion  from  equilibrium 
throat  conditions 

3)  Controls  of  the  integration  to  hit  specific*  area  ratios,  the  nozzle 
throat  point,  the  nozzle  tangent  point,  and  the  requested  end  point 

4)  Controls  of  the  switch  from  the  subsonic  pressure^iefined  equations 
to  the  supersonic  area  defined  equations  when  (M  ^1.02) 

5)  Controls  the  switch  from  implicit  to  explicit  integration. 

For  the  normal  mode  of  operation  of  the  program,  this  subroutine  locates 
the  throat  in  the  following  manner: 

The  gaseous  mass  flow  per  unit  area  (pv)  is  calculated  and  stored  as  a 


: 


* 

* 


i 

4 


function  of  nozzle  axial  location  for  the  present  and  past  integration  step.  When 

(pv)n+l  < (pv)n 

where  n refers  to  the  n*h  integration  step,  the  throat  location  is  calculated  from: 


(X  - x )2- 
n n-i 

(pV)n+i  -(pv)J 

+ (VrV2- 

(pv)  -(pv) 

n n-i. 

2 • 

’(xn+rV 

• r(pv)  -(pv) 

l n 1 

l-l 

"(X  ~X  1> 
n n-1 

•[(Pvkr<r.v)nj 

and  the  n+l**1  Integration  step  is  repeated  using  a step  size  of  X*  - Xr  to  deter- 
mine the  throat  conditions . 

T6  prevent  the  location  of  a false  throat  due  to  roughness  of  an  input 
pressure  table,  ten  integration  steps  are  required  before  the  throat  will  be 


sought. 


Through  the  downstream  throat  radius  of  curvature  the  step  size  is 
controlled  so  as  to  be  less  than  or  equal  to  RWTD*  SIN  (THETA)  /25.0f 
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Subroutine  MAIN1D  also  contains  the  logic  to  control  the  integration  through 
solidification  of  multiple  condensed  liquid  phases  if  they  are  present.  This  logic 
allows  the  integration  procedure  to  hit  the  beginning  and  ending  of  solidification 
exactly  and  to  turn  on  the  solidification  equations  in  subroutine  DERIV  via  the  flag 
IFMELT  in  COMMON  LKMELT. 
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5. 5 >26  Subroutine  OUTPUT 

This  subroutine  provides  conversion  from  internal  computational 
units  to  output  engineering  units  and  the  calculation  of  performance  para- 
meters. The  following  output  parameters  are  computed  by  this  subroutine: 

*« 

The  pressure  (in  PSIA)  is  computed  from: 

P(PSIA)  = p/«33.056 

The  gaseous  species  mole  fractions  are  computed  from: 

R1 

cl.»  ■ T • cl 

The  gas  molecular  weight  is  computed  from: 

Mw  = 49721. 011/lR 

The  percentage  mass  fraction  change  is  computed  from: 

n 

% A (Mass  Fraction)  = 100.0  • (1.0-2  C.) 

i=l  1 

The  gas  heat  capacity  is  computed  from: 

Cp  (BTU/LB-°R)  =3 . 9969  • 10-5  • Cp 
g y 

The  gas  enthalpy  is  computed  from: 

n 

H (BTU/LB)  =3 .9969  • 10-5 ' 2 C,  • h. 

8 1=1  1 1 

At  the  throat,  the  characteristic  exhaust  velocity  (ft/sec)  is  com- 
puted from: 

c*  = A_ 

w p*  • V* 
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Xhe  vacuum  specific  impulse  Is  computed  from: 


ISP. 


V + p/p  - c* 

c _ g- 32.174 


VAC  9 

Ihe  vacuum  thrust  coefficient  Is  computed  from: 


I 


rVAC 


sPvac 

~~C*~ 


Ihe  percentage  enthalpy  change  Is  computed  from. 


%aht  = 


100'  (HREF  -HREF) 

c 


V2/2 


where 


NSP 


HREF  * pl  CA-  ^ + V /2 


HREF  is  HREF  evaluated  at  the  initial  condition  for  the  0DK 
Integration  (i.e.  at  the  initial  contraction  ratio,  ECRAT). 

In  addition  the  actual  molecular  weight  of  the  mixture  of  gases  and  con- 
densed phases  is  also  printed  out. 
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5,5.27  Subroutine  PRNTCK 

For  the  option  to  print  starting  at  step  ND1#  printing  every  ND3rc* 
step  up  to  step  ND2,  this  subroutine  checks  whether  or  not  the  current  step 
should  be  printed.  If  it  is  to  be  printed  this  subroutine  calls  subroutine 
0UTPUT. 
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5.5.28  Subroutine  TAB  GEN  (IFLAG,  L TABLE,  XTAB,  YTAB.  LUSED,  X,  Y,  IERR0R,  N 

This  subroutine  records  a tabular  function  (X,Y(NY))  in  tables  of  fixed 
length.  The  first  and  last  event  will  always  be  tabulated  and  the  table  will 
either  contain  all  of  the  values  specified  or  will  be  at  least  half  full.  Once 
the  number  of  events  exceed  the  table  length,  the  table  will  be  repacked  by 
the  deletion  of  every  other  table  entry  and  tabulation  will  proceed  choosing 
every  2— ^ event (N=0,  1,2...  etc. , where  N is  the  number  of  times  the 
table  must  be  repacked).  The  table  spacing  will  be  a power  of  two  except  for 
the  last  event  which  will  always  be  tabulated. 

The  calling  sequence  parameters  are: 

IFLAG  - denotes  type  of  entry  to  subroutine 
■ - 1 first  entry 
* 0 normal  entry 

*=  + 1 last  entry 

LTABLE  - length  of  tables  available  for  tabulation 
XTAB  - table  for  tabulation  of  the  variable  X 

YTAB  - table  for  tabulation  of  the  variable 

LUSED  - number  of  table  entries  currently  used  (output) 

X the  variable  X 

Y - the  variable  Y 

IERR0R  - error  flag 

NY  - number  of  Y variables  to  be  tabulated 

Note:  One  Dimensional  Mach  Number  Tabulation  Procedure 

At  the  initial  axial  position,  X and  Mach  number  are  recorded. 

TABGEN  is  then  used  with  LTABLE=50  (assuring  25  saved  values).  The  last 
recorded  values  are  the  end  values  for  the  transonic  tables. 
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5,6  Link  50  Subroutines 

This  overlay  and  subroutines  contained  within  are  the  essence  of  the  TD2P 
module.  The  majority  of  this  section  is  from  Reference  4 with  the  major  modifi- 
cations restricted  to  Sections  5,6,2  through  5,6,12.  However,  almost  all  of 
the  subroutines  have  been  modified  to  some  extent  in  order  to  bring  this  calculational 
module  upto  date.  That  is,  this  calculational  module  has  been  updated  from 
F0RTRAN  II  to  F0RTRAN  IV  by  the  extensive  use  of  labeled  common  blocks  and  the 
reductions  in  the  amounts  of  equivalences  to  blank  common  areas. 

Other  modifications  made  to  the  TD2P  module  include  the  extensive  trans- 
mittal of  data  to  reduce  the  amount  and  complexity  of  data  needed  to  execute  this 
module. 


5,6,1  Program  LINK5Q 

This  subprogram  is  mainly  intended  to  facilitate  the  conversion  of  the 
SPP  code  overlay  structure  to  other  computer  systems.  The  only  function  of  this 
routine  is  to  call  subroutine  TD2P, 
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5,6.2  Subroutine  DIST 


Input: 

Output: 


D,  N,  <jg 

u\  * L 

J mT 


M,  ...  N 


Description  : 

Given  an  arbitrarily  large  number  of  spherical  drops  of  mass  mean 
diameter,  D,  which  are  assumed  to  possess  a logarithmic  normal  distribution 
about  D,  then  the  mass  density  is  given  by: 


d D 


[(2  irl/2  D In  <rg] 


exp 


2 


where 

D is  particle  diameter 

m,p  is  the  total  particle  mass 

<7g  is  the  geometric  standard  deviation 


It  follows  that 


The  purpose  of  this  subroutine  is  to  determine  a set  of  N drop  diameters 


j=l,  ...  N 

such  that 


The  method  described  above  is  the  same  as  presented  in  References  10  and  11. 
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5.6.3  Function  EISP 


Function  calculates  the  error  function 

EISP  = c - c (x) 

used  in  finding  the  one  ^mensional  - zero  lag  I as  a function  of  area  ratio,  e. 

sp 

The  independent  variable,  x,  is  either  the  gas  temperature  (before  or  after  solidi- 
fication of  the  condensed  phase)  or  the  fractional  amount  of  solidified  condensed 

phase  (during  solidification).  The  I corresponding  to  the  independent  variable 

sp 

is  also  calculated. 

This  function  is  called  by  subroutine  SEEK  which  is  called  by  subroutine 

ISPID. 
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5,6.4  Function  FAM 


This  function  calculates  the  error  function 
FAM  (M)  = f - C(  M) 
where  c = requested  area  ratio 

€ (M)=  area  ratio  at  the  Mach  number,  M,  given  as  the  argument  of  function 
FAM. 

The  FAM  function  is  called  by  subroutine  SEEK  (Section  5.6.8)  which  is  called 
by  TBLEDG  (Section  5.6.9)  to  determine  the  subsonic-transonic  Mach  numbers 
as  a function  of  area  ratio  for  communication  to  the  TBL  module. 

The  error  function  is  calculated  using  the  following: 

A*/A  = ASA  = *±i  M (1  + xy~  Ms)2^V_1)/(v+1) 

FAM  = e - 1.0/ASA 


5-107 


5,6,5  Subroutine  FIND 


The  purpose  of  this  subroutine  is  to  perform  linear  or  quadratic 
interpolation  from  a table. 

Caning  Sequence 
CALL  FIND  (N,T,X,Y) 
where 

N is  a flag  6uch  that: 

N • 0 implies  linear  interpolation 
N jt  0 implies  quadratic  interpolation 

T is  a table  of  the  format  specified  below 

X is  the  argument 

Y is  the  result 

Usage 

The  table,  T,  must  be  a column  of  a FORTRAN  array.  The  format  of  the  table 
is: 


T(l)  a temporary  cell  used  by  FIND 

must  be  0.  initially. 

T(2)  minimum  argument 

T(  3)  function  of  minimum  argument 


T(2n)  * nth  (maximum)  argument 

T(2n+l)  function  of  the  nth  argument 

T(2n+2)  0.  signifying  the  end  of  table 

T(2n+3)  0.  signifying  the  end  of  table 


If  FIND  is  entered  with  an  argument  above  the  table  maximum  or  below  the  table 
minimum,  it  will  store  the  function  of  the  maximum  or  minimum  argument,  respect- 
ively. If  FIND  is  entered  with  an  argument  which  is  between  the  first  or  last 
two  table  arguments,  linear  interpolation  will  be  used. 
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Functional  Description 


This  subroutine  saves  its  place  In  a table  and>  on  an  option,  either  inter- 
polates by  the  linear  interpolation  formula: 

v . Yg  ' yi  + *.  X,  < * * *0 

*2-*l 

or  performs  a four  point  quadratic  interpolation  by  the  formula: 

2 

y - Afx-Xj)  ♦ ♦ C xQ<xx  < x s Xg  < x^ 

where 

d - *,-*1 

* - V*0 

tmxr*i 

a - (y2  -yj)  (e-*)/d  * y0-gy!  * y3 

d(e-f)  + e2  ♦ ? 

B - (yjj-y^/d  - Ad 

C - yx 

The  above  quadratic  formula  Is  constrained  such  that  the  function  y{x)  passes 
through  points  (x^,  y^)  and  (x,,  yg)  and  has  the  property 

[y(xo)-y0]  - [y(x3)  - y3]  ■ o 


Note:  There  is,  depending  on  the  computer  system  used,  a possible  entry  point 
name  conflict  between  this  routine  and  the  subroutine  of  Section  5.1.6  by  the 
same  name.  These  two  routines  are  not  the  same.  See  Appendix  A for  conversion 

aides. 
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5,6.6  Subroutine  ISPID 


This  subroutine  calculates  ; for  a given  area  ratio  assuming  one  dimen- 
sional equivalent  perfect  gas  flow  with  zero  velocity  and  temperature  lags  between 

the  gas  and  condensed  phases.  This  I is  referred  to  as  I in  most  of  the 

SP  splDO 

documentation  concerning  the  SPP  code.  # 


The  calling  sequence  to  ISPID  is  CALL  ISPID (E, ISP , XI) 

where  E = area  ratio,  c , at  which  The^I  is  desired 

spiDO 

ISP  = I on  exit  from  ISPID  f 

SplDO, 

XI  = initial  guess  at  the  independent  variable 


Before  and  after  solidification  of  the  condensed  phase,  the  gas  temperature 
is  used  as  the  independent  variable  to  iterate  on  the  solution,  however,  during 
solidification,  the  fractional  amount  of  solid  to  total  condensed  phase  is  used. 

This  routine  utilizes  subroutines  SEEK  (Section  5.6.8)  and  EISP  (Section 
5.6.3)  to  determine  the  root  of  the  appropriate  equations. 

The  constants  used  in  these  routines  are  calculated  in  subroutine  AGP 
(Section  5.6.12). 
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5.6.7  Subroutine  SEEK 


This  subroutine  utilizes  subroutine  ITER  to  find  the  root  of  the  equation: 

f 00  = o 


Subroutine  ITER  predicts  new  guesses  of  the  root,  x^,  given  old  values  of  the  error, 
fQ,  and  guess  at  the  root,  xQ,  using  the  method  of  secants.  The  SEEK  subroutine 
' attempts  to  take  maximum  advantage  of  the  properties  of  the  secant  method.  (See 
the  write  up  on  subroutine  ITER). 


Calling  Sequence 


where 

XANS 

X0 

MAXIT 

EPSI 

FCALC 

IB 


B(I) 


Call  SEEK  (XANS,  X0,  MAXIT,  EPSI,  FCALC,  IB,  B,  IT,  IERR) 


the  root  of  f(x)=0  or  the  last  guess  for  x if  the 
maximum  number  of  iterations  has  been  exceeded. 

initial  guess  at  the  value  of  the  root. 

maximum  number  of  iterations  allowed  in  finding 
the  root. 

a so^tion  to  f(x)=0  is  claimed  if 
Jf (x)j  < EPSI  or 

I Ax  | < EPSI/100 . 

an  external  function  subprogram  which  calculates 
the  error,  f(x).  The  calling  sequence  is  F=FCALC(x). 

a bounding  flag  such  that  if 

IB  = 0 no  bounding. 

IB  = 1 the  root  is  assumed  to  lie  within  the 
values  of  B (1 ) and  B(2). 

IB  = -1  the  maximum  absolute  change  allowed 
between  successive  guesses  of  the  root 
is  constrained  to  be  no  larger  than 
|B(l)*|xo|  + B(2)| 

a vector  of  length  2.  The  use  of  this  vector  is  de- 
scribed above. 


IT  = . the  number  of  iterations  required  to  obtain  a solution. 

IERR  = a flag  such  that  if 

IERR=0  a solution  was  obtained. 


IERR=I  failed  to  converge  on  the  root  in  MAXIT 
iterations. 
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5.G.8  Subroutine  TBIEDG 


This  subroutine  calculates  the  subsonic-transonic  edge  conditions  for  the 
trubulent  boundary  layer  module,  TBL.  The  number  of  points,  NTBL,  and  the  initial 
location  of  the  subsonic  start  conditions,  XITBL,  are  determined  by  user  input  in 
the  $TD2  namelist  (Volume  in.  Section  2.9). 


The  method  used  in  constructing  the  edge  tables  for  TBL  is  as  follows: 

1 . The  region  from  XITBL  to  the  MOC  Initial  line  is  broken  up  into  NTBL 
equally  spaced  points. 

2.  Subroutine  SEEK  and  function  FAM  are  used  to  calculate  the  1-D  equi- 
valent perfect  gas  Mach  number,  M _ at  the  area  ratio  where  the 
initial  line  intersects  the  wall, 

K U 

3.  The  axial  position,  Z , in  the  nozzle  which  corresponds  to  an  inlet 
area  ratio  of  nine  is  ftcated. 


4.  The  edge  Mach  numbers  are  then  calculated  from 


M = Min 
e ID 


il  + 


Z-Z 

o 

zpr 

i i o 


M.  -Mini 

■ At  IPlj 
M1DU 


I 


where  the  subscripts 


i i = initial  line 

ID  = 1 dimensional  Mach  numbers  calculated  using  SEEK  and  FAM 
as  a function  of  axial  position,  Z. 


In  addition  to  the  above,  subroutine  TBLEDG  also  recalculates  the  stagna- 
tions pressure  and  temperature  communicated  to  the  TBL  module.  The  values  of 
PQ  and  Tq  are  calculated  from  the  static  P and  T at  the  initial  line  wall  point  in 
order  that  edge  conditions  calculated  by  TBL  would  more  closely  resemble  the 
values  calculated  by  TD2P. 
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5.6.9  Subroutine  TD2P 

Subroutine  TD2P  is  the  master  subroutine  of  the  TD2P  module.  It  controls 
the  initialization  of  variables,  stores  in  the  proper  locations  values  transmitted 
from  other  modules,  reads  in  user  input  in  the  $TD2  namelist,  and  directs  the  flow 
of  the  computation. 

In  the  case  that  particle  sizes  are  not  input  to  the  module,  TD2P  calculates 
the  mean  particle  size  (in  microns)  as 

d = XK  • PCPEXP  • (»)XIEXP  (i  _ EXP  (XL *XL STAR)) 

P 

*(1  + 24  • XD  • r*) 

p = 100  (WPWGT/((1  + WPWGT)  • XMLW)) 

where  XK,  PEXP,  XIEXP,  XL,  and  XD  are  constants  set  in  DATA  statements  or 
they  may  be  input  in  the  $TD2  namelist. 

PC  = chamber  pressure  (psia) 

XL  STAR  = L*  (in) 

£ = moles  of  condensed  phase/lOOgm. 

WPWGT  = m gas/m  condensed  phase 

XMLW  - molecular  weight  of  the  condensed  phase 

The  particle  size  distribution  is  then  calculated  by  calling  subroutine 

DIST. 

Overlay  5,  1 (LINK51)  is  called  to  calculate  the  average  gas  properties 
used  by  the  rest  of  the  TD2P  routines  (see  subroutine  AGP,  Section  5.6.12). 

The  following  items  are  then  computed  and  made  available: 

P 


5-113 


Overlay  5,2  (LINK52)  Is  then  called  for  execution  of  the  Axlsymmetric  Sub- 
sonic-Transonic Subprogram  and  then  the  following  items  are  computed  and  made 
available: 


C 

c 


R) 


N 


P 

r 


m 

P 


(T  R)N 
go 


r* 


r2 

pj 


max 


after  which  Overlay  5,3  (LINK53)  is  called  for  execution  of  the  Axlsymmetric 
Supersonic  Method  of  Characteristics  Subprogram. 


- 


f 


i 
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5,6,10  Program  LINKS  1 


This  subprogram  Is  mainly  intended  to  facilitate  the  conversion  of  the  SPP 
code  overlay  structure  to  the  computer  systems.  The  only  function  of  this  routine 
is  to  call  subroutine  AGP, 
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mm 


5 . 6 . 1 1 Subroutine  AGP 


The  purpose  of  this  subroutine  is  to  calculate  average  gas  properties  based 
on  the  results  of  a one-dimensional  calculation  using  variable  thermodynamic  pro- 
perties (usually , the  results  of  an  0DE  calculation).  These  properties  are  then 

used  in  the  calculation  of  one-dimensional  and  two  dimensional  I in  the  rest 

sp 

of  the  TD2P  module.  The  method  used  is  outlined  below. 


Using  Newton's  method,  an  expansion  coefficient,  v , is  found  such  that 
a polytropic  expansion  from  the  chamber  temperature,  TgQ,  to  the  solidification 
temperature,  T , will  occur  at  the  input  expansion  ratio,  e . The  iteration  equa- 

pm  *nm  ^ 


tions  used  are 


where 


pm 


— (l+l)  _ — (i)  /c  cpmx  /d 

V ~ V l 


'<f  - » 
V 


1/2 


+ 1 


dy 


QO 


pm, 


_X+J_ 

2(y  - 1) 


and 


(?-]) 


It 

d y 

Next  the  value  of  the  specific  heat,  C , for  the  gas  phase  is  calculated 
from  the  one  dimensional  energy  equation  as 


U 4 


qm 


pg 


2V  V 


W \ 

w 

A 

P 

. j 

• 

wg/ 

Wg 

'pi 


i.e.  the  reference  enthalpy  at  stagnation  is  taken  as  zero.  The  flow  velocity 

at  the  T is  given.  It  is  also  assumed  that  C . for  T ^ T _ is  also  known, 
pm  3 pi  pm 

An  expansion  coefficient  for  the  gas  phase  is  next  calculated  assuming 


no  lag: 


so  that 


1 - (y  - 1)  cPi 


iat  C 
9 pg 


Next  the  value  of  C^,  the  particle  specific  heat  after  solidifi- 
cation (T  < Tpm)#  is  calcuiated  from  the  energy-equation  as 


w (1*  - T ) 

p'  pm  Aps' 


i 1 + 


u - c fr  - t ) 

gs  wpg  ' go  ps' 


Cpl  “ V - -- 

wg 


The  values  of  T and  U are  input  values  corresponding  to  some  point 
ps  gs 

after  solidification.  The  above  value  of  Cps  is  then  checked  against  the  actual 

value  of  C at  the  solidification  temperature.  If  the  calculated  value  is  more 
ps 

than  10%  different  than  the  actual  value,  the  actual  value  is  used.  The  reasons 
for  the  above  test  are 

a)  it  prevents  ridiculous  values  of  C from  being  used 

ps 

b)  it  prevents  the  results  of  the  TD2?  calculation  from  being  too  sensi- 

tive to  the  area  ratio  selected  at  which  values  of  T „ and  U „ are  ob- 
tained ps  gs 

If  the  Prandtl  number,  Pr,  is  not  input  or  transmitted  from  the  0DE  module, 
it  is  calculated,  using  kinetic  theory,  as 


Pr  = 4V/(9V  - 5) 
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5.6.12  Program  LINK52 


This  subprogram  is  mainly  intended  to  facilitate  the  conversion  of  the  SPP 
code  overlay  structure  to  other  computer  systems.  The  only  function  of  this  rou- 
tine is  to  call  subroutine  PARTIL. 
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5,6.13  Subroutine  ABCALC 


Coefficients  are  calculated  for  first  and  second  order  transonic 
approximations : 


First  Order; 
C - 


u 2 - 1 
0 


- k-1  2 

1 ~ TTT  U 
k+1  O 

C a 
o 


20  AC. 


u a 
o 


21  2C, 


V21 0 


40  8C, 


Second  Order; 


A201 

=5 

, k-l 

2i5ruo a 

4A201  a20  + Co  S 

b21 

VC 

4C1 

A102 

B 

2 

a + u 3 
0 

* 

2uo  a6  + A102  a + 

4A201  a21 

b22 

m 

4C1 

A120 

b 

2u  b01 
0 21 

A220 

0 k-l 

2 k+T  uo  a21 

A320 

k+r  u0  a2o 

2Co  b22  + A120  a 

+ 4A220  a20  + 2A320  a21 

b40 

16CX 

A121 

=3 

2a21  “ + 4uo  b22 

*321  “ k+1  uo  a21 

2uoa218  + 4uob22°  + A121a  + 4A220a21  + 16A201a40  + 2A321a21 
b41  “ 16C1 

A140  “ a212+2uob41 
A340  “ k+f  Uo  a40 


jA 


5.6414  Subroutine  CCALC 

Coefficients  are  calculated  for  the  third  order  transonic  approxi- 
mations: 


Third  Order; 


B 

k-1  A 

°202 

k+1  A102 

c 

Co  y/2  + 4A201  b21 

+ 4B202  a20 

C22 

4C1 

B103 

■ u i + aB 
0 j 

c 

4A201  b22  + 4B202  3 

21  + A1026  + B103a  + uoay 

C23 

4C, 

n , k-1  2 

B120  " 4 k+1  a20 


t.  k-1 

*220  " k+1  A120 


*40 


B12Qa  + 4B220  a20  * 2A320  b21  * 2Co  C22 
16C, 


B121  ” 4uo  C22  + 2b2l“  + 8 k+1  a20  a21 


r k-1  , 

°221  k+1  121 


B321  ■ iST  (a20“  + uo  b21> 


'41 


16C.  (4A320  b22  + 2A321  b21  + 16A201  b40  + 4A220  b21 


+ B121a  + 4B221  a2Q  + 4B22Q  + 2B321  a21  + A^B 


+ 6Co  C23  + 4uo  C22c) 


"122 


6uo  C23  + a216  + 4b22°  + 4 M a2l' 


(a01a  + u b00) 


*42 


B 


"80 


16CX  (B122°  * 4B221a21  " 16B202a40  * 2B322a2*  * A1216 


+ 2A102b22  + 16A201b41  * 4A220b22  * 4A321b2? 


12uoc23a  + «oan  y) 


k-l 


140 


2uoC41  + 2a21b21  + 16  k+I  a20a40 


k W . 

a240  * k+1  A140 


B340  * k+1  (s20  a21  + 2uo  b40} 


C60  " 36C1  (B140a  + 4B240a20  15B220a40  * 2B340a21  + 2A12Qb22 


+ 16A220bAo  + 2A340b21  * 4A320b4i,  + 2CoC42  + 4uoa21C22* 


k-l 


141 


4a21b22  + 2b41a  + 16  k+1  a21a40  * 4uoC42 


B34l  " k+1  (a21  + 2a40°  + 2uo  D41) 


JL  , 


+ ^A^'>i^ai  + ^6A*n  b-n  * * 


“61  36CX  v "34CU22  321^41  ^w201  60  >20  43 


+ fc14la  + 4B240a2i  + l6E221a40  + 2B341a21  * A140£ 


2A121b22  + 4uoC42a  + l2uza21C23> 


160 


2a  C,..  + 2a01  b. , + 16  T7~  a,n2 
o 61  21  41  k+1 


B360  • k?T  (2a21  a40  + 3uo  b60} 


64C1  v 16.T 


+ A6B240a4C  + 23360a21  + 2A140b22  + 36A>™b‘ 


220  60 


~ 4A340D41  * 4uoa21C42) 
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5,6.15  Subroutine  DCALC 


A dummy  subroutine  containing  only  a RETURN  statement.  This  sub- 
routine is  reserved  for  the  calculation  of  coefficients  for  the  fourth 
order  transonic  approximations. 
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5.6.16  Subroutine  FCALC 


This  subroutine  calculates  the  boundary  condition  functions,  f 
from  coefficients  and  derivatives  determined  by  subroutines  ABCALC, 
CCALC,  and  PCALC: 


u - u 


t dr 

f0  * u — - v 
2 o dx 

2 

f du  dr  , d r _ dv 

*3  dx  dx  Uo  , 2 ” dx 

dx 

d2u  dr  - du  d2r 

t . * — - — + 2 r + u - 

4 . 2 dx  dx  , 2 o,3 

dx  dx  dx 


,3  ,2 

dr  d v 


dx4 


.3  j ,2  ,2  , .3  ,4  ,3 

d u dr,  - d u d r . du  d r d r d v 

, 3 dx  + J . 2 , 2 J dx  , 3 o , 4 ~ 3 

dx . dx  ax  dx  dx  dx 


If  a special  throat  expansion  (see  Step  3,  subroutine  PARTIL)  is 


being  performed,  the  following  variables  are  computed: 


* , - o , o 

u ■ u +a  x +8  — + yt~ 
o o 2 '6 

2 

x 


<**  - a + fix  +Y72 
o L 


If  transonic  conditions  in  the  nozzle  inlet  are  to  be  fared  to  the 

throat  conditions,  the  following  variables  are  computed: 

2 3 

. , . (z  +x  ) . (z  +x  ) 

A k fJk  y n k 'wo 

U *=  U +•  a (z  +x  ) + £ - + Y 7 

c wo  2 6 


* O*  * (zw+X0} 

a +3  (z  +x  ) + Y ^ 

WO  2 
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5,6,17  Subroutine  JAMES 

This  subroutine  provides  control  for  the  perfect  gas  transonic  flow 
approximation  computations.  Newton's  Method  is  employed  to  find  solutions 
(s<*e  subroutine  NEWT) 


, sonic  point  displacement 
UQ  > axial  velocity/critical  flow  velocity 

a , 1st  order  coefficient 

6 , 2nd  " 

Y , 3rd  " 

to  the  non-linear  equations  (see  subroutine  FCALC) 


Solutions  are  sought  such  that: 

if  1st  order  then  3,  y,  f^  = o,  and  f^  » o are  deleted 
if  2nd  order  then  y and  fg  = o are  deleted 

if  the  conical  inlet  angle,  q.,  Is  greater  than  the  faring  angle,  then  a 
special  expansion  is  performed  at  the  nozzle  throat  with  uQ  and  T*,=o  deleted. 
Expansions  at  points  (rw,zw)=l  upstream  of  the  faring  point  are  performed 
with  uQ  and computed  by  subroutine  FCALC  and  fg  set  zero. 


Calling  Sequence 

RR  rw,  radial  wall  position  coordinate 

XX  zw,  axial  wall  position  coordinate 

TT  Vrc»  rec^Proca^  throat  radius  of  curvature 

VEK  k,  the  average  gas-particle  expansion  coefficient 
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5.6,18  Subroutine  LEGS 


Purpose 


TO  solve  an  over  de  termined  linear  system  of  equations  AX  * b in  the  least 

T T 

squares  sense,  i.e.,  to  solve  the  "normal  equations"  A AX  » A b. 

The  superscript  T denotes  the  transposed  matrix. 


Usage 

Calling  Sequence: 

CALL  LEGS(AB,X,BKDS,M,N,MD,ND,ATA,R,Z,SUS, SUSP, II, 12,13, 14,I5,IDLE) 


where 

AB  * the  M by  N variably  dimensioned  augmented  matrix  (A,b)  with 
maximum  row  dimension  MD  and  maximum  column  dimension  ND. 

X * the  solution  vector  of  dimension  N. 

BNDS  * the  bounds  vector  of  dimension  N. 

M = the  number  of  rows  in  the  augmented  matrix  (A,b). 

N - the  number  of  columns  in  matrix  A. 

MD  m the  row  dimension  of  array  AB.  MD  i M. 

ND  * the  column  dimension  of  array  AB.  ND  £ N + 1. 

ATA  = a one  dimensional  array  in  which  the  upper  triangular  part  of  the 
augmented  normal  matrix 

ATA  ATb  \ 

bTA  bTb  J 

is  stored  by  rows.  The  dimension  of  this  array  mutt  be  greater 
or  equal  to  (N+l)  (N+"’)/2. 

R * a one  dimensional  array  used  for  storage.  If  the  inverse  of  aTA 
is  requested,  the  lower  triangular  part  is  .tored  by  ~ows  in 
this  array.  The  dimension  must  be  greater  than  cr  equal  to 
((N+2)  (N+3)/2)  - 1. 

Z c an  output  indicator.  If  the  bounding  -“lion  is  losen,  L t 
indicates  that  the  solution  was  affected  by  the  be  nds. 

SUS  ■ the  squared  norm  of  b,  i.*,,  SUS  « ||b||  ■ 


5 - 1 ? 7 


SUSP 


p 

■ the  squared  norm  of  AX-b,  i.e.,  SUSP  ■ | |AX-d| | 

U « Bounding  option.  If  II  • 0 the  solution  is  unbounded  and  the 
routine  minimizes  | |AX-b  1 1 . Fbr  II  / 0,  the  routine  minimizes 
||AX-b||2  under  the  side  condition  that 

2(Xi/Bi)  * 1# 

The  sum  extends  over  all  i for  which  BA  > 0.  If  BjL  • 0,  the 

routine  ignores  the  i'th  row  and  column  of  the  augmented  normal 
matrix  and  sets  Xi  « 0.  If  Bi  < 0,  the  routine  does  not  bound 

XA»  The  constants  B ^ must  be  stored  in  ENDS.  When  II  « 0, 

the  ENDS  vector  is  set  to  -1. 

12  m Solution  option.  If  12  ■ 0,  the  routine  solves  for  X.  If 

12  £ 0,  the  BNDS  vector  is  placed  into  the  X vector  and  SUS, 

T »1 

SUSP,  and  (A  A)  (if  requested)  are  computed.  This  may  be 
used  for  checkout  purposes  when  the  routine  is  used  in  an  iterative 
process  for  solving  nonlinear  systems.  .It  is  also  used  for 
comparing  two  "solutions"  of  a nearly  singular  system. 

13  « ATA  option.  If  13  « 0,  the  augmented  norm  A matrix  is  computed 

and  placed  in  ATA.  If  13  > 0,  the  normal  matrix  computation 
will  be  omitted.  It  is  assumed  the  user  has  already  calculated 
this  matrix  and  placed  it  in  ATA.  If  13  < 0,  the  normal  matrix 
is  computed  and  added  to  whatever  is  in  ATA.  This  is  helpful 
in  case  the  augmented  matrix  must  be  partitioned. 

Ik  m Inverse  option.  If  14  / 0,  the  normal  matrix  is  not  inverted. 

Upon  exit  R will  contain 

a.  the  lower  triangular  part  of  the  matrix  S,  stored  by  rows, 
where  F is  an  Nth  order  lower  triangular  matrix  with  (-1) 
on  the  main  diagonal 

b.  the  solution  vector  X 

c.  -1 

d.  the  main  diagonal  of  the  diagonal  matrix  D. 

Let  C m A^A  + zB  ^ , then  S and  D have  the  property  that  SCS^*  « D* 

If  Ik  « 0,  the  matrix  C is  inverted  as  C"1  * STD‘1S  and  C"1 
replaces  S before  exit  from  the  routine. 

15  * ATA  only  option.  If  15  / 0,  the  routine  exits  after  forming  the 
normal  matrix.  If  15  * 0,  the  routine  continues  by  checking  12. 


* A more  complete  explanation  is  found  under  Functional  Description. 
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IDLE  * ATA  row  deletion  option.  In  normal  use  IDLE  Is  a single  cell 
. integer  of  value  0 and  has  no  effect  on  the  solution.  The  use 
of  this  option  is  best  explained  by  an  example.  Suppose  the 
user  forms  the  augmented  matrix  and  wishes  to  delete  certain 
rows  from  consideration  in  the  least  squares  problems  (because 
of  bad  data  or  some  other  reason).  He  wishes  to  delete  rows 
20#  31-35#  5 and  l4~l6.  This  can  be  accomplished  by  putting  the 
following  sequence  of  integer  into  IDLE:  20,0,31#35#5#0,lU,l6,0 

The  list  is  terminated  by  a zero  in  an  odd  position.  IDLE  must 
be  dimensioned  in  the  driving  program  if  it  is  to  be  used. 


Functional  Description 

1.  Bounds  - When  the  bounds  option  is  chosen  the  program  minim*  the 
quadratic  form  ||AXrb||2  subject  to  the  side  condition  £ 1, 

-2  -2 

where  B is  the  diagonal  matrix  whose  i’th  diagonal  element  is  B. 

if  B^  > 0,  whereas  if  B^^  < 0,  the  i'th  diagonal  element  is  equal  to 
zero  (if  B^  ■ 0,  the  routine  ignores  the  i'th  row  and  column  of  the 
augmented  normal  matrix  and  sets  **  0). 

The  method  of  Lagrange  multipliers  is  applied  to  solve  this  problem 
Thus,  the  system  of  .linear  equations 

(ATA  + zB*2)  X = ATb 

is  solved  for  various  choices  of  the  real  positive  multiplier  z,  until 
a solution  X( z ) is  found  which  satisfies  the  side  conditions. 

The  search  for  an  appropriate  value  of  z is  systematized  as  follows: 

With  6 * l/lO  and  £g  " 2/l0,  the  routine 

o 

a.  Finds  X ■ X(0).  If  (B~  X,  X)  £ 1 +6,#  the  problem  is  considered 
solved.  Otherwise, 

b.  define  y(z)  ■ (B*^  X(z),  X(z)).  Now,  y(o)  > 1 ♦ €^.  Compute 
y(z,  ),y(z2),...,  vhere  z^  ■ z ^ + 10(j-l)H  and  H « M|X 

[B^(  |ATb)i|/B;.  - ATA(i,i))]  whenever  this  expression  i6  positive, 
otherwise  H « M^X  [.0001  ATA(i,i)B^] 

The  computation  of  the  y(z,)*s  is  continued  until  a z is  found 
such  that  1 «tg  £ y(zj)  S1  + ^ in  which  case  X(z^)  is  accepted 
as  the  solution,  or  until  two  values  of  z^,  renamed  z^  and  Zg 
are  found  such  that  y(z^)  >1  + and  y(zg)  > 1 - £g. 

The  required  value  of  z is  now  bracketed.  Now 

c.  a value  of  z^  is  chosen  with  Zg  < z^  < z^.  If  1 «£g  £ y(z^)  £ 1 + £. 
X(z^)  is  accepted  as  the  solution.  Otherwise, 
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inverse  quadratic  interpolation  to  zero  is  applied  to  obtain  a 
new  estimate  z^.  If  1 - < y( z^)  S 1 + X(zu)  is  accepted 

as  the  solution.  Otherwise, 

e.  select  from  the  set  and  zh  that  pair  which  bracket  the 

solution  most  tightly  and  return  to  (c). 

The  routine  will  halt  this  iterative  process  if  the  number  of  solutions 
of  the  linear  system  reaches  twenty. 

T -2 

2.  Solution  of  the  linear  system  - Let  C ■ A A + zB  . The  routine 

T 

finds  two  matrices  S and  D such  that  SCS  ■ D,  where  D is  diagonal 
and  S lower  triangular  with  minus  ones  on  the  main  diagonal.  It 
is  easy  to  find  S and  D for  a 1 x 1 matrix  C.  Suppose  they  hs.ve 
been  found  for  a k x k matrix  C.  Now  augment  C by  another  row  and 
column 


and  seek  a vector  CO  and  a scalar  3 such  that 

o\  /c  d^  /sT  w \ /d  o\ 

„t  ij  (d »*;  (o  j ■ (o  t)  ■ 


It  16  easy  to  verify  thatuli 

B- 


STD'1Sd 
2 -M>Td 


6a+isfy  the  requirements.  The  routine  builds  the  matrices  S and  D 
by  the  above  process  the  final  result  being  a decomposition  of  the 
augmented  matrix. 

/ S 0 \ /aTA+z3"2  ATb  W ST  \ Id  0 ) 

-lj  \ bTA  bTb/  yo  -1/  (o  c</' 

The  vector^? is  the  solution  vector.  The  above  method  is  a variant 
of  the  "square  root"  method. 


3*  Residual  Calculation  - The  residual  is  calculated  as  follows: 


||AX  - b||2  - (ATAX,X)  - A(ATb,X)  + (b,b)  . 
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5.6.19  Subroutine  NEWT 


Given  a set  of  m functions  in  n unknowns  (m  * n) 
f1(x1,  x2,  ...»  xn)  1 - 1,  ...  m 

this  routine  will  attempt  by  Newton’s  method  to  find  values 

Xl’  X2*  Xn 

such  that 


® 2 
l X 

1-1  1 


Is  a minimum. 


Restrictions 

The  user  must  supply  initial  values  for  the  solution  vector  and  also 
a subroutine  to  evaluate  the  functions.  The  subroutine  must  communicate 
with  NEWT  through  C0MM0N  statements.  Subroutine  TRW  LEGS  is  required. 

No  error  exits  are  provided. 

Method 

Newton’s  method  is  used  to  iterate  for  a solution  vector.  The  matrix 
inversion  is  performed  by  subroutine  TRW  LEGS. 

In  the  event  Newton’s  method  yields  a vector  which  is  farther  from  a 
solution  in  a least  squares  sense  than  the  previous  estimate,  the  increment 
vector  is  halved. 


Usage 

Calling  Sequence: 

CALL  NEWT (M , F , MF , N , VAR , NVAR , FCALC , WF , PERTV , EPS1 ,KAXIT , 
NUMIT,BIN,AB,B,KB) 

Index  to  Calling  Sequence: 

Input:  M, MF,N, VAR, NVAR, WF, PERTV, EPS1 ,MAXIT, B, KB 

Output:  VAR,F,NUMIT 

Working:  BIN,AB 
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Explanation  of  Calling  Sequence: 


1) 

2) 

3) 


4) 

5) 

6) 


7) 

8) 


9) 


10) 


ID 

12) 

13) 


M is  the  total  number  of  functions,  m. 

F is  the  array  of  the  m function  values  found  by  FCALC. 

MF  is  an  array  of  m control  words: 

if  MF(I)«0,  include  F(I) 
if  MF(I)«1,  exclude  F(I) 

N is  the  total  number  of  variables,  n. 

VAR  is  an  array  of  the  n variables,  x^  j * l,...n 

NVAR  is  an  array  of  n control  words: 

if  NVAR(J)*0,  include  VAR(J) 
if  NVAR(J)*1,  exclude  VAR(J) 

FCALC  is  a name  for  the  subroutine  which  calculates  the 
functions,  F(I).  A program  which  calls  NEWT  must  have  an 
EXTERNAL  statement  containing  this  name. 

WF  is  an  array  of  m weighting  factors.  These  are  used  in 
conjunction  with  EPS1  to  determine  whether  a solution  has 
been  reached. 

a>1  - WF(I) 


PERTV  is  a perturbation  factor  used  in  calculating  the 
partial  derivatives  required  by  Newton's  method.  If 
PERTV  « s,  then: 


3F 

3x 


F (x  4 6)-  F (x) 
6 


_2 

where  5 = max( | ex | ,e*10  ; 


EPS1  is  ac  error  bound  used  to  determine  whether  a 
solution  has  been  reached.  If  = EPS1,  then  a solution 
is  claimed  when 


m 

Z 

i*l 


W 


MAXIT  is  the  maximum  number  of  iterations  allowed. 

NUMIT  is  the  number  of  iterations  required  for  solution. 

BIN  is  an  array  used  internally  by  the  routine  and  must  be 
at  least  ot  dimension  n^  + 7n  + 3. 
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14)  AB  is  a tvo-dimensional  array  containing  the  augmented  matrix 
of  subroutine  LEGS  and  must  be  at  least  of  dimension  AB (m,  n+1). 

15)  B is  an  input  bounds  vector  of  dimension  n.  If  the  bounds 
option  is  chosen,  the  iteration  increments  Ax  are  constrained 
to  be  less  than  B^  (see  LEGS).  B is  not  modified  internally. 

16)  KB  is  an  input  flag  to  be  set  equal  to  1 if  the  bounds  option 
is  desired.  Otherwise,  set  KB  - 0.  (See  LEGS.) 
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5.6.20  Subroutine  flNED 

The  relations  for  one-dimensional  two-phase  perfect  gas  nozzle  flow 
•re  integrated  through  the  nozzle  conical  inlet  (Region  I of  Figure  5 6-2) 

The  coordinate  system  is  centered  at  the  apex  of  the  conical  inlet.  The 
differential  equations  are: 


dV 

-L 

dz' 


f ' r*  (V  - v ) 
8 p1 

2 m * = '~ 


P r. 


J * •**  j 


max 


dT 

dz* 


- 3 ,s  g;  r. 


. r2 

P P 


J 


P cn 
r P, 


(T  - T ) 

_2j L 

V 

pj 


j"l»  •••j 


max 


dV 

£. 

dz 1 


V y 



V - yRT 
8 g 


r2RT 


+RCpiJmax  VdTPj 

Cp  fa1  ; ST 

8 g 


where* : 
P 


j w 
Jmax  n 

£ 

J"1  w 8 

8 8 


dV  n 

' V V 


o m UJV  A 
g g g 


* 2*r*z  (1  - cos  0i)2 


,2 


T ■ T 

g 2C 


V - jmax  wp 


-r-  L ^ 


Pg  J’1  L * PJ  > 


CP  - T ) + 


■] 


" Y 


'T  iN 


V 


* The  area,  A is  taken  as  a spherical  cap  of  half  angle  0 
and  radius  r*z'.  6 i 
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1+2.52  C 


f* 

P 


1 


J 


1+3.78  C 


1/ 


PJ  PJ  l 


1 + 


3.42C  f 
J P4 


y1/2  P 
g * 


f (R  ) (by  interpolation  from  Table  5 • 6—*l) 
P 


2p  r (V  - V )/y 
g Pj  8 P£ 

y /p  r /~RT 


/ 


g g P 


j 


g 


The  initial  conditions  are: 


200 


199+Y  8o 


1 + W - x> ! 


^max 

i+  E — 

■I*1 s_ 


W 


C 1 

P 'max  p 

1+r^S  r 

V*1  8 


CD  j w 
P Jmax  p 

1 + tl  + (Y.  - 1)  ^ B]  Y -r-1 


C 


and 


zo  ■ - ir<Ao/2,T(1  ' cos  9i)>1/2 


A 

o 


m 


200 


199-hr 


i/(7-D 


/ p 
go  g0 


The  integration  proceeds  until: 


z'  - - (1  + R (1  - cos  6.))/sin  e. 

& C il 


using  an  integration  step  size 


h - Az  , which  is  the  input  item  DZMIN 
min 


The  one-dimensional  gas-particle  velocity  lags, 
are  determined  at  position  z* 


and  thermal  lags,  L 


V 


The  quantities 


and 

u / u*  = Vg/u* 

s 

are  also  calculated  to  be  used  by  subroutine  FCALC, 
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5.6.21  Subroutine  PARTH 


This  subroutine  is  the  master  subprogram  for  the  Subsonic-Trjnsonic  Link. 
Tne  axisymmetric  gas-particle  flow  field  of  Region  II,  Figure  5.6-1  is  computed 
using  the  procedure  outlined  below. 

STEP  1: 


The  variables  and  are  computed; 

1/2 


equil 


2CP  Tg 
- 8 O 


/ kequil  -1  \ 
kequil  +1  ' 


z * -R  sin  0 . -R  ( 1-cos  0.) 

SCI  1 


R * [ 1 + Rfi  (1-cos  0j)  ] / sin  0^ 


initially 


k 


equil 


STEP  2 

Subroutine  0NE  D is  used  to  determine  conditions  at  the  upstream 

end  (z  = z ) of  Region  II,  Figure  5.6-1 
s 


Region  II  is  next  divided  into  zones  as  illustrated  by  the  example 
given  in  Figure  5.6-2  JWall  coordinates  for  the  perturbation  points  are: 

r * 1 + R (1-cos  0_J 
w c P 

z “ R sin  0_ 
w c P 


where 


V1 


0 


j 


— 0 . 9 **  ~ 0J»  • .0,  , **.  0, 

i nu  i nj  j 


i.e.  a total  of  n.  + n,  +1  zones  are  constructed, 
i 3 


STEP  3: 


The  perfect  gas  transonic  flow  subroutine,  JAMES,  is  used  to 
determine  the  variables  xo,  uq,  a,  8,  y and  the  transonic  approximation 

coefficients  for  each  zone. 


If  faring  is  necessary  (0^  > 0^)  a special  throat  expansion  is  first 
performed  and  the  variables  u*,  a*,  y*  and  0*  are  computed. 

STEP  4: 

One-dimensional  gas-particle  flow  relations  are  integrated  along 
the  nozzle  axis  from  zg  to  zax^g«  The  differential  equations  integrated 

are  those  given  in  subroutine  TRACE.  Initial  values  for  the  gas  variables 

T /T  and  u are  evaluated  using  subroutine  PR0P  Initial  values  for  the 
g g0  g 

particle  variables  are  computed  from  and  as  determined  by  subroutine 
0NE  D.  Transonic  zones  are  exchanged  as  soon  as 


, UQi  * U°i  ± 1 

u * 2 

g 

Values  K , and  L.  at  z . are  made  available  for  use  in  Step  8. 
j j axis  r 

STEP  5: 


A corrected  estimate  (constant  lag)  for  the  expansion  coefficient, 
k,  is  computed  using  values  at  z » zaxis 


k « 1 + 


JCdL 


1 + c2y 


•^max  pj 
°2  " i i ~ 

j - 1 g 


w r c 


LJ  + KJ  (1"V  C1 

L g. 


CP  jmax  V 

i + r1  Z er1  LJ 

1 w 
Jmax  woj 

1+  Z — Kj 

j-1  g 3 
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the  variable  u * is  recomputed  as 


Using  the  above  expansion  coefficient,  k,  transonic  conditions  near  the 

initial  line  are  recomputed  by  subroutine  JAMES.  A corrected  estimate 

for  the  nozzle  mass  flow,  m , is  found  by  integrating  gas  properties 

S 

along  the  initial  supersonic  data  line  using  subroutine  WDGI. 


STEP  6: 


Steps  2 through  5 are  repeated  twice  to  assure  convergence  in  the 

variables  k and  m . 

g 


STEP  7: 

ei  ei  ei 

For  9n  " 20  ’ 3~  * J~  • and  6i 
an  initial  position  of 


r R sin  0 
o n 


z - z + R (1-cos  0 ) 
os  n 


R * [l  + Rc  (1-cos  0i)  ] f sin  0^ 

is  selected  on  the  upstream  boundary  of  Region  II,  Figure  5,6-2,  and 
particle  trajectories  are  integrated  through  Region  II  using  subroutine 
TRACE. 


Tables  are  constructed  for  the  following  variables  as  found  along 
the  initial  line: 


J »n 


® (1-cos  0 ) 

Pj  n 

2 Tir*^  (1-cos  0^) 


K 


!s 


j.n  u 


g 
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K 

j.n 


j,n 

2 

Cj,n 


g 

T - T 
8o  Pj 

"T  -t 


where  n - 1 , . . 4 
j - 1...J 


max 


Limiting  particle  streamline  trajectories  correspond  to  n **  4. 


STEP  8: 

An  initial  supersonic  data  line  is  constructed  consisting  of 
NILP  + 1 points  with  equal  radial  spacing.  This  data  line  extends  from 
the  axis  to  the  wall  along  the  parabola 


axis 


_/faxis_l Jwall  j r2 
rwall 


Limiting  particle  streamline  points  are  added  along  this  data  line  and 
those  equally  spaced  points  which  fall  too  close  ^within  '^NILP  ) 


are  deleted. 


Gas  properties  P , p , u , and  v are  computed  at  each  point  using 
6 g g g 


subroutine  PR0P. 


Particle  properties  u , v , h , and  p are  interpolated  at  each 

pj  pj  pj  pj 


point  from  the  tables  constructed  in  Step  4 and  Step  7.  The  interpolation 
formulas  are: 
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where 


P(x  )-x  + a(x.  ) r2+b(x  )r4+c(x  )r6 

j.n  j.o  J»“  j.n  o,n 


c(x  )-  — r 
J»n  iJ1 


<rj,4  ‘ rJ,3)<XJ,2  - Xj,o> 


,2  x 2,2  2 w 2 2 N 

(rJ,«  “ rj »3  rj,2(rj,3  " rJ,2  <rj,4  " 'j^ 


(x  - - x ) 
j.3  j,o 


(XJ.4  " Xj,o> 


2 2 2 2 2 2 
rJ*3  (rJ,3  " rj,2^  rj»4  <*J,4  - rJ,2> 


b(x.  ) - — 5 5 

J »n  (r2  _ rZ  ) 

J.3  rj,2; 


Xj .2  ~ x1.o 
2 


c(x.  )(r?  ~ + r2  ,) 

J.n  J.2  j,J 


a(x.  „) 
J.n 


XJ,2  " X3,o 


- b(x  )r  . 

J.n  j.Z 


c(x,  )r  , 
J.n  j.2 


The  subscript  n*o  denotes  axis  values.  The  variables  K.  and 

J t° 

LJ  are  made  available  in  Step  4 while  K*  , and  \l>.  are  taken  as: 
j»o  r j.o  j.o 

K'j,o-K'j.l 

♦j.o-’j.l 

STEP  9: 


For  each  initial  supersonic  data  line  point  the  varables  p , u , 

8 8 


v , r,  z,  h , p , u , and  v are  printed. 

8 pj  pj  pj  pj 
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'‘I  : ■ . - . .* 


5.6.2  2 Subroutine  PCALC 


The  following  quantities  are  computed: 


x ■ x , r * r , x * z 
o w w w 


dr  ,22 

S " xw  <VV> 


d^r  2 2 ? 2 ° “3^ 

— r - (R  -x  ) + x2  (R  2-x  *) 

^x2  c w wcw 


d3r  2 2 “3/2  3 2 2 ‘5/2 

“ 3x  (R  -x  ) + 3x  J (R  -x  l) 

d^3  w c w w c w 


d^r  2 2 ~3^2  2 2 2 4 2 2 

» 3 (R  -x  ; + 18x„  (R  -x  Z)  + 15x4  (R  2-x2) 

,4  cw  wcw  wcw 


dx 


u*  «ffx  + r* 


u - 6 ■=■  +(b01+2b00x)  r + b,,  r 


2 21  22 
3 


41 


2n  2 


Y 7 + (2c22x+3c23x  ) r + (cyi+2c/0x)  r + cci  r 


41  42 


“61 


u + uf  + u*  * + u* ' 1 
o 


2(a20+a2lx)  r + 4a40 


2(b21x+b22x2)  r + 4 (b40+b41x)  ^ + 6b60  r* 


2(c22x2+c23x3)  r + 4 (c40+c41x+c42x2)  r3 


+ 6 (c60+c61x)  r5  + 8c80  r? 


v'  + v"  + v'" 
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and 


du 

dx 


, dr 
u + u "7“ 
x r dx 


^ - u + 2u  — + u + u (j1)2  + 3u  Og£-) 
.2  xx  xr  dx  r , 2 rr  \dx/  rrx  ax 
dx  dx 


2 

dx2 


/dr  \2 


/dr  *3  . Qll  dr  d r 
unr  (3F)  + 3urr  cGT 


dx 


A 

dx3 


u + 3u  ~ + 3u  — — 4 u — r 

xxx  xxr  dx  xr  , 2 r . 3 

dx  dx 


dv 

dx 


. dr 

V + V “7” 

x r dx 


d2v 

dx2 


A 

dx3 


V + 2v  £ + v + v f£| 

xx  xr  dx  r dx2  rr  \ dx / 


v + 3v  — + 3v  + v — ~ + 3v 

xxx  xxr  dx  xr  dx2  r dx3  rrx 

+ v {l1)  3 + 3v  LJi 

rrr  \dx|  rr  d'  dx2 


I dr  \2 

\ dx  / 


where  the  partial  derivatives  u , u , u ,u  , u ,u  %u  . \ 

x*  r xx*  xr  rr’  xxx*  xxr* 

u , v , v . v , v , V vV  .v  .v  , and  v have  been 

rrr’  x*  r*  xx*  xr*  rr'  xxx*  xxr*  rrx*  rrr 

by  differentiation  of  the  expressions  for  u and  v given  above. 


rrx’ 

obtained 
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5.6.23  Subroutine  PRgjP 


Given  the  coordinate  position  (r,z)  the  gas  properties  u , v , 

T / T , and  p are  calculated  from  transonic  coefficients  previously 

g g0  g 

determined.  The  velocity  derivitives  u , v , u , and  v , may  also 

8z  gz  8r  8r 

be  computed. 

The  quantity  I in  the  calling  sequence  is  a flag  such  that: 

1*1  implies  u , v T /T  , and  p are  to  be  calculated 
g g g gQ  g 

1=2  implies  u , v , u , and  v , are  to  be  calculated 

O O O <3 


gz  gr 


k-1 
’3  “ k+1 


q.  * z_  + x - Z * z 
V 1 o w 

2 3 

P = u + a z + B + y T~ 
o o 2 6 


r = o,  Properties 


u * u P 
g g o 

v =0 

f M „*2 


g„ 


g 


P - Pe  =»  t-1 

g 6°  Tgo 


r f o Properties 

* a21  + ^21  + ^^22  + C22^  z + z 

P4  = b41  + CA1  + 2(b42  + C42)  ’ f 


P,  * b,  + c.  + 2c.  0z  + 3c.  z 
6 61  61  62  t>3 


P8  * C81  f 2c82Z 
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P10  " C101 

P1  " a20  + (a21  + b21)  2 + (b22  + c22)  ^ + c23z3 


a40  + b 40  + c40  + (b41  + c41)  z + (b^  * c">>  z + c,,z 


42  42' 


"43 


b*n  + c*n  cti)  z "*■  (c__)  z "*■  c£oz 

2 


5 “60  ' “60  ' '“61  ’ “61'  ~ ‘ "-62'  ‘ ' w63 


7 b80  + c80  + C81Z  + c82z 


P9  “ c100  + c101z 

P11  “ c120 
Q4-r2 


* r 2 4 6 8 in  l 

Ug“Ug  LFo  + P2r  +V  + V +P8-  +P10r  ] 

vg  - 2 ug*  [ Pxr  + 2 P3r3  + 3 P5r5  + U P?r7  + 5 P?r9  + 6 Pnr 

2 2 
k-1  , u + v 

g 


ii 


K-l  / U TV  v 


P„  “ 

g g, 


a ) " 


r a o.  derivatives 


G ■ a + Bz  + y — 
O L 


u ■ u G 

g2  g o 

v • 2 u * P, 

gr  g 1 


u - 0 


v ■ 0 

g. 
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r t o,  derivatives 

3 

2 


u - u * f 2 ?,  r + 4 P,.  r3  + 6 P^  r5  + 8 P„  r7  + 10  P,„  r9 
8^.  S L 


10 


v *=  u 


r2  + 30  Pc  r4  + 56  P,  r6  + 90  P„  r8 


\-"8*[2Pl*12P3-  •"■S'  •-•7 

+ 132  P1Qr10  1 


G2  ■ 2(b22  + c22)  + 6 C23Z 


G4  = 2(b42  + C423  + 6 C43Z 


G6  " 2 c62  + 6 C63  2 


G8  = 2 c82 


6Z  g 


u [ G + G,  r2  + C,  r4  + G,  r6  + GQr8  1 
g I o i 4 t 8 
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5,6.24  Subroutine  TRACE  

This  subroutine  integrates  particle  trajectories  through  Region  II 
of  Figure  5.6-2  using  the  differential  equations: 


where 


drJ 

V 

Pi 

dz 

du 

pj 

<“s 

-v 

dz 

- “ a — 
0 

dv 

pj 

V 

dz 

- ■ a 

O u 

pj 

dT 

pj. 

-b 

0 

V 

dz 

cpt 

u 

2 V’  /J 

.N  f 

A _L 

a * 
0 

2 a \ T 

P l 

;o  r2p 

L _ 

a 2 CP* 

D * 
0 

ao  3 p 

r 

'■'I 

with  f'  and  g'  as  defined  in  subroutine  0NED. 

1 J 
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Initial  conditions  are  computed  by  subroutine  PARTIL  using  subroutine 
0NED.  The  initial  position  (r^  * r^,  zq)  is  selected  on  the  upstream 

boundary  of  Region  I,  Figure  5.6-1. 

The  integration  is  terminated  with  the  intersection  of  the  particle 
trajectories  and  the  end  line  of  Region  II,  Figure 5. *6-1. This  end  line 
is  the  parabola 


/Zaxis  Zwall  \ 2 

'axis  “l  2 ) r- 


wall 
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5.6.25  Subroutine  WPG! 


This  subroutine  integrates  along  the  initial  supersonic  data  l'ne  for 
mass  flow  (units  of  lbs  mass)  by  trapezoidal  rule: 


* “n  *2 
*5+1  - 4 


zi+l  ‘ zi 

P u + p U - ■ (p  v +0  V 1 

8i+l  8i+l  8i  8i  ri+!  ~ ri  81+1  g1+1  gt  g± 


The  Initial  data  line  is  a parabola  defined  as 


with 


ri+i  “ rt  - -02 


z - z -I  Za]tls  ~ zwall  12 

i+1  axis  l 2 Jri+1 

‘wall 


rl  " 1 


i - 1,2,. ...50 


Particle  mass  flows  are  also  computed: 


I w 


m 


\ * 


max 
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5.6.26  Program  LINKS 3 


This  subprogram  is  mainly  intended  to  facilitate  the  conversion  of  the  SPP 
code  overlay  structure  to  other  computer  systems.  The  only  function  of  this  routine 
is  to  call  subroutine  N2MAIN. 
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5.6,27  Subroutine  ACOMP 


Given  two  adjacent  points,  m and  n,  whose  properties  are  known,  the 
Reynold's  number  associated  with  the  jth  particle  size  at  a third  point  in 
the  neighborhood  is  assumed  to  be: 


where 


A drag  coefficient  is  found  by  quadratic  (or  linear)  interpolation  on 
Reynold's  number  from  Table  5.6-1. 


f - f (Re  ). 
P P 3 


A slipflow  correction  is  applied  to  the  drag  coefficient; 

2 


(1+9. 45c)  (l+2.52c)+3.03t/ 


PJ  p 1(1+9. 45c)  (1+3. 78c)+.9U(4+11.34c)c2 


1 + 


(r#j 


C f 


where 


(u  -u  + (v  -v 

8 Pj 8 Pj 


Re, 


RT 


The  coefficient,  A , is  defined  as 

pj  * 

r 

O 8- 


f' 


PJ  2 (T  R)N 
8o 


a. 


. r2 
P P, 
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5«6«28  Subroutine  ADJK 


This  subroutine  performs  Step  4 as  described  in  subroutine  KPBPT.  The 
two  arguments  in  its  calling  sequence,  CALL  ADJK(BARL , Q)  are 


3ARL  « X input 

Q ■ output 
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5 >6. 29  Subroutine  AXISPT 

This  sub.?utlne  uses  two  known  points  to  calculate  a new  axis  point 


Properties  at  point  3 are  computed  from  properties  at  points  1 and  5. 
The  gas  and  particle  streamlines  passing  through  point  3 lie  along  the  axis. 

The  following  steps  specify  an  axis  point  calculation.  Tests  are  made 
to  avoid  particle  calculations  when  no  particles  are  present. 

STEP  1 and  STEP  2 

Each  variable  at  point  3 is  set  equal  to  the  average  of  the  correspond- 
ing variable  at  points  1 and  3 except: 


STEP  3 

Previous  values  of  point  3 are  saved,  then: 
X-  T A ( Ar  X3) 
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STEP  4 


! 


U “ “ + “T*  (*3  " z5)(Bi  + Bi  > 

PJ3  PJ5  2 3 5 J3  J5 


K -K  + 2 <23  • 25> (E1  + E1  > 

Pj3  PJ5  2 3 5 J3  J5 

2p  v 

P1  P1 
J5  J5 


Pj  u 

J3  Pj. 


P u 

P1  P1 
. J5  J5 


<Zo  * 


*3> 


STEP  5 


2<J3  - Jx) 

V L1  -S-G3 


Vg; 

Pg3  ’ J3  " 2 


P ■ P + 

83  g5 


(S3+S5)(Pg3-  Pg5)-(T3+T5)(23-z5) 


Steps  3 through  5 are  iterated  until  the  calculated  variables  change 
by  less  than  a prescribed  amount,  or  until  a maximum  number  of  iterations 
is  reached. 
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5.6.30  Subroutine  CHECK 


This  subroutine  compares  current  and  previous  values  of  the  point  3 


flow  variables  P.p.u.vh  o ,,  „ . 

8 8 8 8 Pj*  PP j’  P d Vp  for  relat*ve  con- 

vergence by  calling  subroutine  GRIT.  If  . variable's  encountered  failing 
the  test  the  program  returns  with  the  convergence  flag  set  as  failed 

(N0  - 1).  If  convergence  is  achieved  the  subroutine  returns  with  the 
convergence  flag  set  as  passed  (N0  - 0). 
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5-6,31  Subroutine  CNTRL 

The  purpose  of  this  subroutine  Is  to  control  the  construction  of  the 
finite  difference  mesh  for  the  method  of  characteristics  solution  of  the 
supersonic  nozzle  flow.  Left  running  characteristics  are  constructed 
starting  at  the  nozzle  throat  point  (r  ■ 1,  z - 0).  These  left  running 
characteristics  extend  from  the  Initial  data  line  or  nozzle  axis  to  the 
wall*  The  mesh  points  are  calculated  by  subroutines  INPT,  AXISPT,  WLPT 
and  KPBPT  under  control  of  this  subroutine. 

Additional  points  are  inserted  In  the  mesh  by  subroutine  CNTRL  by 
property  averaging.  Points  are  inserted  to  control  the  maximum  mesh 
spacing.  In  a nozzle  of  shallow  radius  of  curvature  at  the  throat,  the 
limiting  particle  streamlines  may  be  very  closely  spaced.  A point 
Insertion  and  drop  technique  allows  this  situation  to  be  handled  without 
the  introduction  of  an  unnecessarily  large  number  of  character 1st ice. 

No  provisions  have  been  made  for  the  crossing  of  particle  trajectories  or 
for  the  occurrence  of  shocks.  Weak  shocks,  however,  will  be  ignored  by 
point  deletion  or  characteristic  mapping.  Figure5*  6-3fllustrates  a typical 
mesh  construction  containing  Kth  particle  boundary  points,  interior  points. 
Initial  line  points,  and  points  inserted  and  deleted. 

Whenever  a mesh  point  calculation  or  insertion  is  completed  success- 
fully, the  output  routine,  PRINT,  is  called  so  that  the  point  may  be  printed. 

Subroutine  CNTRL  also  integrates  the  wall  pressure  by  trapezoidal  rule 
to  determine  the  axial  component  of  force  existing  between  adjacent  wall 
points, 

r2 

4F  - 2nr*2  / P„  rdr  - £(P  + P )(r,2  - r 2)r*2  . 

r5  8 *■  g2  g5  1 5 

Total  thrust  is  then  found  by 
2 12 

F.  ■ F + 2irr*  / P rdr 

1 i 8 

where  F is  the  thrust  across  the  initial  line  as  calculated  by  subroutine 
C0NSTS. 
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The  quantities  defined  below  are  not  computed  in  subroutine  CNTRL 
but  are  basic  to  the  point  calculations  controlled  by  subroutine  CNTRL. 
These  quantities  are  given  here  for  convenience. 


TP. 


- 


/p  (u  2 + v 2)  - Y? 
g g g g 


u p - yP 
g g g 

p u v + .R  -R 

&._B— 8 L-2- 


P*  U«  V*  ~ 1R  2R 
1V 
V 

-R 

u 

g 


P<  u 


P4 

j 

p 'jN  (u  - u ) 

!&  _E Ii_ 

P u 

■gj  pi 

f (v  - v ) 

>!*!  1 !i. 


1} 


j 


it 

P8J 


yP  V 
t,t  - — 


VVg  1R  + U,  2R) 

1R  xv 


°*(V»  1»  ~ U*  2R) 
XR  XV 
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Pg  ♦ <F3  w3  + P2  M2)  + ^3R2  + 3*3> 


u 


(6  8 

{ ‘ ~r  <pg2  vg2  + pg3  vg3> + ~r  \ \ + pg3  ug3> 


(r3  - r2) 


(*3  " z2) 


p 

N 

P 

g3 

f + 

g2 

P , 

2 3 

P 

g3 

g2 

N 


2*2J 


P 

83 

N 

P 

g2 

iPJ»  i 

83 

3*3  + 

k i 

g2 

N 


3t2l 


+ 4 t«g2(r3  - r2)  - y^z3  - «2) 


% '2  * *2> 


8-5 


'p  ^ 
8 


i.p»  J 

g3 


1*3  + 


lV  } 


J,  - P 

3 8i 


<z3  " 21> 


2 (f3  n3  + f3  n3)  - jtjRj  + 3n3) 


V + Pft  v ) + -r—  (p  u + p u ) 
81  g3  83  2 81  81  83  83 


(r3  - V 


(Z3  - 21) 


p 

g3 

N 

P 

81 

lPo  J 
g3 

2*3  + 

k i 

gi 

P 

g3 

N 

f 

P 

81 

ipE  J 

g3 

3*3  + 

k i 

gl 

2*ll 


N 


3tl] 


+ T Iu8l(r3  - V - vgl(z3  - 21>  - ug3  rl  + vg3  *il 


fp  1 

8’ 

N 

,t.  + 

N 

Pgl 

11 

V 
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Note:  When  calculating  an  axis  point  the  evaluation  of  reoulres  an 

Indeterminate  (v  - r * 0)  quantity,  . In  this  case  the 
S ^ 

expression  below  is  used: 


► « 

V 

u 

V 

i 

00 

g3 

r 

* a 

3 rl 

V 

81 

L u + (z,  - z,)  - — J 

v 

This  calculation  amounts  to  extrapolating  for  the  ratio,  assuming 
that  the  flow  near  the  axis  is  a source  flow. 
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5.6.32  Subroutine  CONSTS 


An  Integration  Is  performed  using  trapezoidal  rule  for  thrust 
along  the  initial  line: 


F = 2„r*2 


j-  r 

-wall  ( -v  ) + T,  0 u (u  -v  ) rdr, 

J [ g Hg  g g r j i pPj  Pj  Pj  Pj 


where  K is  the  total  number  of  particle  sizes  present  at  the  point  under 
consideration. 

The  indeterminate  quantity 

v 

r 

at  the  initial  line  axis  point,  m,  is  calculated  as 


& 


lx) 

r j 


^m-1 

Vl 


g, 


m 


L Vl 


+ (z  - z ) 
m m-1  r 


*m-l 


m-1  J 


This  quantity  is  required  in  both  the  axis  and  off  axis  point  cal- 
culations. 
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5,6.33  Subroutine  CRIT 


The  purpose  of  this  subroutine  is  to  compare  tvo  values  for  relative 
convergence  and  return  a flag  stating  if  convergence  has  been  achieved. 

Calling  Sequence: 


XN 

m 

NO 


m 


A flag  such  that 


if  x - 0 then  NO  - 0 
n 

if  |(x_  - x )/x  |s  5.10’5  then  NO  * 0 
m n n 

if  |(x  - x )/x  |>  5.10“5  then  NO  « 1 

m n n 
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5.6.34  Subroutine  CUBIC  (X.  Y.  VP,  N,  ARG,  YARG) 

The  purpose  of  this  subroutine  is  to  perform  cubic  interpolation 
for  a tabulate  1 function  whose  derivatives  are  known. 

Calling  Sequence; 

X 

Y 

YP 

.N 

ARG 
YARG 

Restrictions: 

The  calling  program  must  define  arrays  for  the  dummy  variables 
X,  Y#  and  YP.  These  arrays  must  be  at  least  of  lengths  N + i9  N#  and 
N respectively.  The  subroutine  will  save  its  last  used  table  position 
number  in  X(N  + 1). 

If  x < the  program  returns  y = y^ 

If  x > Xjj  the  program  returns  y = y^ 

Mfilhsd: 

Given:  x ^x<x, 
o 1 

so  that  yQ#  yo'#  y^,  and  y^1  are  known.  The  cubic  interpolation  formula 
given  below  is  used  to  determine  y 

y = A(x  - + B(x  • XQ)^  + C(x  - xq)  + D 

where 

a = ^ ((y  i'  + y0‘) h - 2kI 

B-prlfri+2^.,)h-3kI 


independent  variable9  x^9  such  that 

dependent  variable9  y.,  = y(x^} 
derivitives  of  the  dependent  variable9 

Yi  s 

is  the  number  of  entries  in  each  of  the  above  tables; 
is  i#...  N 

is  the  argufttent,  x9  for  which  interpolation  is  requested 
is  the  result9  y = y(x) 


is  a table  of  the 

*i+l  ixi 

is  a table  of  the 

is  a table  of  the 
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D s y 

• Jo 

kaU-r0 
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5.6-35  Function  EFN 


The  quantity  is  calculated  In  this  function  subroutine.  A test  is 
made  on  particle  enthalpy,  h , in  order  to  determine  the  proper  particle 


phase. 


v 


u 


N 


u 


N 


K * V1  ■ " 

I1-.  - r%] 


•bp  R)  g1  r 
pi  I p / * 

J \ at 


h - h p 

Is Ll  _ A. 

Cp  R p 


h > h 

pi  p/ 


h < h < h 
ps  Pj 


h < h 


Pj  Ps 
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5,6,38  Subroutine  ERR0R 


The  purpose  of  this  subroutine  is  to  print  error  messages  as  requested 
by  subroutines  COTRL,  ADJK,  and  WLPT.  See  Appendix  A for  conversion  aides. 
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5.6.37  Subroutine  KP3PT 


This  subroutine  uses  five  known  points  to  calculate  a Kth  particle 
boundary. 


Properties  at  point  3 are  computed  from  properties  at  point  2 and 
interpolated  properties  at  point  1. 

The  following  steps  specify  a Kth  particle  boundary  point  calculation. 
Tests  are  made  to  simplify  the  calculations  when  a gas  and  one  r article 
boundary  point  is  to  be  computed. 

Step  1 

Each  variable  at  point  3 is  set  equal  to  the  average  of  the  corres- 
ponding variable  at  points  2 and  4.  Point  1 is  initially  set  equal  to 
point  4. 
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Then 


J ■ 1|  •••* 
j - 1,  K-l 


Step  2 

Previous  vslues  of  point  3 are  saved,  then 


7 - TA(<2>  *3) 


TA(C  , C ) 
\ K3 


r3  “ r2  "*■  k^z3  " z2^ 


Step  3 

I - ta(x3>  xx) 

Points  6 and  7 are  selected  such  that  point  6 is  point  4 and  point  7 
is  the  next  point  on  the  previous  characteristic. 
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Step  4 


This  step  is  performed  by  subroutine  ADJK. 


If  qx  > 1 
If  < 0 
If  1 > qx  > 0 


select  new  points  6 and  7 which  are  one  point 
higher  on  the  previous  left  running  characteristic 
and  restart  Step  4. 

select  new  points  6 and  7 which  are  one  point  lower 
on  the  previous  left  running  characteristic  and 
restart  Step  4. 
proceed . 


Step  5 


r6  + 1l(r7  ’ r6} 


P + q.  (P  - P ) 

86  1 87  86 

Po  + - P»  > 

g6  1 g7  86 

u + q-  (u  - u ) 

86  1 87  g6 


v + q,  (v  - v ) 
g6  1 g7  g6 


h + q,(h  - h ) 
P<  1 P<  Pj 


V * - \ > 

J6  37  36 


j ■ 1 K-l 

J - 1 K-l 


u + q,  (u  - u ) 
P4  P4  P4 


j * 1,  ....  K-l 


v - v + q,  (>r  - v ) 

% \ 1 % \ 


j * 1*  •••»  K-l 


Step  6 

Steps  6 and  7 are  performed  for  each  j - 1,  . . . , K.  The  gas  prop- 
erties at  point  4 are  first  saved. 


If  J - K 


If  j * K 


restore  gas  properties  at  point  4 and  go  to 

Step  7. 

proceed. 


TA(C  , C ) 

p4  p4 

h h 


r3  + V(z2  " z3  - r2  T, 


1 - C 


: Zl~Zz 
Pj  - r2 


r4  ~ r2 
rl  ' r2 


Z4  “ z2  + VZ1  - z2> 


Pg  + q4(Pg  " Pe  > 
g2  81  g2 


■ \ *vv  V 

■ \ * - V 

• \ * v«4l  - y 


V * "‘‘V  - V > 

J2  31  J2 
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V • pP,  + Vv  - pp.  > 


PJA  PJ 


Jl  ri. 


Step  7 


Step  8 
for  j - K 


u - u + q.  (u  - u ) 

'j*  »j,  * "j,  -j. 


v ■ v 4*  q.  (v  v ) 

P1  pi  4 P1 

J4  32  J2 


' 1 M 

J; 


% vh 


K + + B,  ) 


2-3  ' J j 

+ 2L(z3  - V(Dj3  + Dj4) 


hPj  +2<Z3"  V(EJ3  + V 


Pj3  lV  <r32-  r22)  " vp  (23-z3)(r2 


rj- 


rJ. 


I P [u  (r2-  r 2) 

+ r7)l  ' % PJ4  4 2 


- V ^-z2)(r2  + VI  +0n  tu  (rA2-r32) 

J4  J2  J2 

- V I(Z4  - z2)(r4  + r2>  - (z3  - z2)(r2  + r3)1]  } 

J o 


for  J < K 


V (tl2  - r22) 

J3 


J — { 

- V [U3  - Z2)(r2  + r3)  + (zx  - z,)^  + r3)J  { 

J 1 


Pn  [u„  (r 

P1  P1 
J4  J4 


(r3  - ) - v [(z4  - z2)(r4  + + 0^  - z4>(r4  + r )]] 
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1 “ z4)(r4  + rl)  “ (*1  " z3)(rl  + r3)H 
z4  - z2)(rA  + r2)  - <z3  - z2)(r3  + r2)]]} 


j “ 1 K-l 

J - 1 K-l 

J “ 1 X-l 


.1  - 1 K-l 


A test  Is  cade  to  determine  whether  point  5 Is  above  or  below  the 
Xth  particle  streamline. 

If 


then 

h - p m u *=  0 and  primes  must  be  added  in 

% % *5 


Step 

10 

to  Jv  G3,  Gj,  Hj,  Hj,  T3,  and  Tjj 

otherwise: 

h 

y(h 

+ h ) 

% 

2 PK 

4 

PK 

*2 

P 

a. 

f(p 

+ P ) 

% 

2 PK 
*4 

PK 

*2 

u 

m 

“(n 

+ u ) 

PKc 

2 PK, 

PK0 

5 

4 

2 

V 

m 

|(v 

+ v ) 

% 

2 PK 
*4 

PK 

K2 

Step  10 

For  the  primed  quantities  the  implied  summations  are  to  be  taken  for 
j ■ 1,  . . . , K-l,  not  K. 

2(Jj  - J2MQ3  + Q2  + h5  + a3)  - 2(J1  - j2hq*  + q3  + Qj  + q2) 

Ug3  “ (L’+  L3+  L|+  L2)(Q3+  Q;+  H5+  H3)-(L3+  L2+  G5+  G3)(Q’+  Q3+  Qj+  Q2> 


2(j;  - J2)  - ug  (Lj  + L3  + + L2) 

Qj  + Q3  + + Q2 


J2  + 


Ug3(L3  + V + + V 


2 


<S3  + 


Sq 


- P 


S5)(Pg3 


" V " 


(T, 


+ T5)(*3  - 


V 


8« 


Steps  2 through  10  are  Iterated  until  the  calculated  variables  change 
by  less  than  a prescribed  amount  or  until  a maximum  number  of  Iterations  is 
reached. 


5,6,38  Subroutine  N2MAIN 


The  purpose  of  this  subroutine  is  to  call  the  two  portions  of  the  axisymmetric 
supersonic  method  of  characteristics  subprogram  (subroutines  N3MAIN  and  CNTRL) . 


5,6,39  Subroutine  N3MAIN 

The  purpose  of  this  subroutine  is  to  call  subroutines  C0NSTS  and  WALL. 
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5.6,40  Subroutine  PRINT 


The  purpose  of  this  subroutine  Is  to  print  output  for  the  supersonic 
flow  calculation  as  requested  by  the  program  input.  Printout  may  occur 
after  the  completion  of  each  mesh  point  calculation.  Points  for  print  are 
selected  as  follows: 

1)  The  following  points  are  always  printed: 

axis  points 

Kth  particle  boundary  poluts 
initial  line  points 
wall  points. 

2)  Interior  points  are  selected  for  print  only  along  every  n^th 
left  running  characteristic  and  only  at  every  n^th  position 
along  these  characteristics. 

3)  Inserted  points  are  printed  if  all  points  are  to  be 
printed  (n^  * n0  ■ 1). 

The  items  printed  are  listed  in  the  table  below  in  the  order  they 
appear,  left  to  right,  on  the  output  sheet.  A header  is  printed  for 
identification  purposes  above  each  characteristic. 

Row  One: 


Item 

Header 

Meaning 

Units 

L.R.C.  number 

LRC 

left  running  characteristic  number 

none 

Xdent.  number 

ID 

type  of  point  (see  below) 

none 

r 

R 

r position  coordinate 

none 

z 

Z 

z position  coordinate 

none 

M 

MACH 

Mach  number 

none 

T 

g 

TG 

gas  temperature 

°R 

V 

ft 

VG 

gas  velocity  (scalar) 

ft/sec 
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Row  One:  (Cont'd) 


Item 

Header 

Meaning 

Units 

6 

g 

THETA-G 

streamline  angle 

degrees 

T /T 
8 *0 

TG/TGO 

ratio  of  gas  temperature  to  chamber 
temperature 

none 

p /p 
8 80 

PG/PGO 

ratio  of  gas  pressure  to  chamber 
pressure 

none 

P„/P0 
8 8o 

DG/DGO 

ratio  of  gas  density  to  chamber 
density 

none 

V*  P /P 
L ?J  8 

SDK/DG 

ratio  of  total  particle  density  to 
gas  density 

none 

CF 

CF 

thrust  coefficient 

none 

I 

sp 

ISP 

specific  impulse 

sec. 

iteration  10. 

IT 

number  of  iterations  required 

none 

Rows  Two  through  K+l 

A row  is  printed  for  each  particle  size,  k>l,...,K. 


Item 

Header 

Meaning 

Units 

k 

K 

particle  size  number 

none 

**k 

REK 

particle  Reynolds  number 

none 

V 

Pi, 

VPK 

particle  velocity  (scalar) 

ft/sec 

K 

e 

Pk 

THETA-K 

particle  streamline  angle 

degrees 

K 

T 

Pk 

P„  /p„ 
Pk  8 

TPK 

particle  temperature 

°R 

DPK/DG 

ratio  of  particle  density  to  gas 
density 

none 

Pn  /0r, 

pk  po 

DPK/DPO 

ratio  of  particle  density  to  chamber 
particle  density 

none 

r 

Pv 

RPK 

particle  radius 

ft. 

■i 
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Of  the  quantities  above,  the  following  are  calculated  within  sub- 


routine PRINT: 


T, " pg  R 


V * * u + v 

g g g 


e - ate  tan  | • 57.29578 

8 \V 


2 P r . /p  T R 
« Pk  8 80 

■■k*  ”57- 

80  \ 8 


(u  - u )2  + (v  - V )2 

8 Pfc  8 pk 


'*  *2 

"t  p. 


SP  w + 8 v 

8 E Pi 


h - h 

pk  p/ 
T - T + ~r L 

pk  pm  CP / 


for  h > h 
pk  p £ 


T - T 
pk  pm 


, for  h ^ h ^ h 
ps  pk  p/ 


ft,  C, 


, for  h < h 

pk  ps 
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The  table  below  relates  the  printed  Identification  number  to  the  type 


of  point  calculated: 

ID 

1 

2 

3 

4 

5 

6 
7 


Type  of  Point  Calculated 

initial  line  point 
interior  point 
axis  point 

Kth  particle  boundary  point 
wall  point 

point  inserted  on  the  previous  L.R.C, 
point  inserted  on  the  R.R*C. 


In  addition,  whenever  a wall  point  is  printed,  subroutine ISP1D  is  called  to 
calculate  the  ID  reference  ISP  and  the  quantity,  is  calculated  as 


I (r  ,Z) 

spTD2P  w 
*TD2P  “ I (r^J 

sPlDO A w 


5.6.41  Subroutine  PUNT 


This  subroutine  uses  three  known  points  to  calculate  an  Interior 
point. 


Properties  at  point  3 are  computed  from  properties  at  points  1, 
2,  and  4.  Point  4'  is  required  only  for  the  calculation  of  particle 
density,  Ppj^»  **  point  4'  is  not  available  (IFG  ■ 0)  an  alternate 
calculation  is  used  (see  STEP  6). 


The  following  steps  specify  an  interior  point  calculation.  Tests 
are  made  to  avoid  particle  calculations  when  no  particles  are  present. 

STEP  1 

Each  variable  at  point  3 is  set  equal  to  the  average  of  the 
corresponding  variable  at  points  1 and  2.  Other  values  required  for  the 
first  iterations  are  set  as: 


3 

3 


j-1 K 
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STEP  2 


Previous  values  of  point  3 are  saved,  then 
k • T A ( <2$ 

7 - T A ( Ax,  X3) 


rl  - r2  + < *2  “ * zi 


< - X 


STEP  3 


tj  - t2  + « u3  - i;> 


crIA  <v  v 


r3  * °g  [ 


*r*2 


z2  " z3  ' r2  r1-r2 


r5  ‘ 


Vz2 


1-  C 


8 Vr2 


r5-r2 


rl_r2 


Z5  “ Z2  + q5  (zrz2> 

P - P 4 q.  (P  - P ) 
g5  82  5 8X  82 


Pg5  “ P82  + q5  < P8l  - Pg2) 
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\ ■ \ * 'i\  ■ v 


Vo  ’ V„  + 1s<v.  " ▼_  ) 

85  82  ^ 8^  82 


hPj  “V  + q5(hp  - V } J“l- 

J5  J2  h J2 


Pn  “ P„  + *c<P « “ P„  ) 

p4  P4  5 p p, 

3 5 ^2  J2  J2 


V-V*VV-V>  w- 


15  J2 


II  J2 


V ■ V + Vv  -V  > i-1*” 

35  J2  J1  j2 


Steps  A,  5,  and  6 are  repeated  for  each 


C - T A (C  , C ) 
P4  P1  P1 

J J3  h 


r„  + C 


r _ zi  " z2  1 

3 Pj[  2 3 2 rx  - r2  J 


- zi  - z2 

1 - C — 1 


PJ  rl  " r2 


u-x-^ 

A r — r 
rl  r2 


*4  ’ *2  + V21  " z2> 


P„  - P + q. (P  - P ) 

84  82  * 8X  g2 
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\ - X * "*'%  - V 


u • U + fl.  (u  - u ) 
«4  *2  4 *1  *2 


V “V  +Q/(v  - V ) 

*4  *2  4 *1  g2 


hp4  “ V + q4(hP.  ■ hP,  ) iml K 

34  32  31  32 


pPl  " pp,  + vv  • pPl > J-1 K 

34  J2  31  J2 


u -u  + q.(u  “ u ) j-l,...,K 

P*  pi  4 P1  pi 

34  J2  31  J2 


v -v  + q.  (v  -v  ) 1-1,..., K 

P,  P.  HAV  p p 

h J2  h 


Calculate 


V - V+— <*3- V <V  +Bp  > 

J3  J4  j3  J4 


v - v + -r (a.  - z,)(D  + D ) 

P1  P1  2 3 A Pi 

J3  J4  J3  34 


V “V  +2(*3-V(V  + Ep  > 


’3  J4 


J3  "34 
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STEP  6 


If  point  4 ' is  not  available  (IFG  • 0) : 


. <ri2  - 4>  - v [iuL  + (ti + r 

33  33  L S 


3)(zrY 


Y Y (rl  ‘ r3)_  Y I (Z1  ' l2)(rl  + r2> 
J2  J2  J2  L 


2 2 
r3  " r2 


/ 2 2v 

Y r3  -r2> 


-vp  (zrz2)(r1+r2)  - (r1+r2)(zrz3) 
Ji  L 


If  point  4 * is  available  (IFG>0): 


[ <V*2>  (r2+r3)  + (zrZ3)(rl+r3)] 

Jo 


Y (rrr2)-vPj  [ (V'Z2)(V  +r2)  + <zr*4,)(r4,+rl)" 


> u <rj-r*)-v 

Pj4.  PJ4.  PJ4. 


L 

+Y  Y (r3-r4,)_vp  [ (V,4')(r*'+rl)-<*l-*3><rl+r3)] 

J1  J1  31 

fY(r4’'r3)'Y  [ (,«'-*2)(r4'+r2)-<*3_*2)(r3+r2)] 

J2  J2  J2 


2(J3-J2)  (Q3-H)2+H5H3)-2(J1-J2)  (2Q3-H31-K)2) 

53  tfLj+L^)  (Q3+<32+H5+H3)-(L3+L2+G5+G?)  (2Q3+Q1+Q2) 
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2(J3-J2)-ug  (2L3+Lx+L2) 


V83  “ 2Q3^2 


P„  " P„  + 
*3  85 


“g3 

2 

(W(?,i-F,<)-(iJ«5)(,3.,5> 

2 


Step  2 through  Step  7 are  iterated  until  the  calculated  variables 
change  by  less  than  a prescribed  amount,  or  until  a maximum  number  of 
iterations  is  reached. 
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5.6.42  Subroutine  SUMP1 


Routine  SUMP1  calculates  the  qualities  it  where 


1*  » T)  A P 

j-l  PJ  PJ 


2‘  - E v v u. 


j=i  pj  pj  ?j 


f'  A P v 

j=i  PJ  PJ  PJ 


k r j 2 

(y-i)  z An  p„  (U  ■ -U  ) + (v  -V  ) - 

J-1  pj  PJ  L 8 vi  J 


C E u 
p 

3 Pr 


u 1 
i-il 

(if 

iW 


5.6.43  Subroutine  SUMP2 

Subroutine  SUMP2  calculates  the  same  quantities  as  described  above  in 
SUMP1 . Routine  SUMP2  also  calculates  the  additional  quantities  l*  where  K is 
replaced  by  K-l,  which  are  required  by  the  particle  boundary  point  calculation. 
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5-6.44  Subroutine  TAFN 


This  subroutine  calculates  an  angle  jveraged  tangent  using  the 
relation: 


r* 


TA  ■ tan  I *r  (arc  tan  X + arc  tan 


Y)] 
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5,6,45  Subroutine  WALL 


This  subroutine  makes  available  five  options  for  the  specification  of  nozzle 
wall  contours.  Either  a cone,  parabola,  arc,  contour,  or,  a spline  fit  curve  may 
be  attached  tangentially  to  a circular  arc.  The  circular  arc  is  specified  as  follows: 

Symbol  Definition 

Angle  defining  the  downstream  end  of  the  arc 
0j  Angle  defining  the  upstream  end  of  the  arc 

r/r*  Normalized  arc  radius 

where 


One  hundred  points  are  generated  along  the  arc  for  equally  incremented 
values  of  angle,  q: 

r = 1 + r/r*  (1  - cos  q)  , 0L  = 9*  0j 
x = r/r*  sin  ^ 

Five  options  are  specified  below  for  attaching  wall  contours  tangetially 
to  this  circular  arc.  The  cone  option  requires  the  nozzle  expansion  ratio  as  In- 
put. The  parabola  and  arc  options  require  the  exit  point  coordinates,  r and  x , 

v C 

as  input.  The  spline  fit  option  requires  a smooth  table  of  wall  coordinates,  r{ 
and  x.,  and  the  exit  angle  at  the  nozzle  lip,  o , as  input. 

j e 
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Option  1,  Cone 


xc  = (re  - ra)ctn  0 + xg 

Option  2,  Parabola 

Two  hundred  points  are  generated  at  equally  incremented  values  of  x such 

that: 

r = a + v b(x  - c) 


where 


*%  - 1 ~ 2ra(xe  ~ xa}  tan 
2 Cre  ' ra  ' Ue  " V tan  »j] 


b = 2(r  - a)  tan 


c = x - 
e 


(r  - a)2 
e 


Option  3,  Arc 

Two  hundred  points  are  generated  at  equally  incrmented  values  of  angle 
such  that: 

r = ra  + R~(cos  q - cos  e^) 
x = xa  + R(sin  0^  - sin  q) 

where 

R = 


(x  - x )2  + (x  - r )2 
e a e a 


T[(xe  - xa)  sin  6j  - (re  - ra)  cos  ^ 
and  the  increment  on  angle  is 


200 


where 


»e  = Sin 


■*[- 


sin  9j  - 


X (xe  - xaJ 


'] 


Option  4,  Spline  Fit 

Two  hundred  points  will  be  placed  along  a curve  input  in  tabular  ft  mi  by 
a cubic  spline  fit  method  (see  subroutines  CUBIC  and  SLP) . Derivatives  are  cal- 
culated from  the  input  values  for  0.  and  the  exit  nozzle  lip  angle,  0 . This  curve 

J 6 

must  connect  tangentially  to  the  circular  arc  described  above. 
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5.6.46  Subroutine  WLPT 


This  subroutine  uses  two  known  points  to  calculate  a wall  point. 


Properties  at  point  3 are  computed  from  properties  at  points  2 and 
5.  The  wall  is  necessarily  a gas  streamline.  It  is  assumed  that  no 
particles  are  present  at  points  2 and  5. 

STEP  1 

Each  variable  at  point  3 is  set  equal  to  the  average  of  the 
corresponding  variable  at  points  2 and  5. 

STEP  2 

Previous  values  of  point  3 are  saved,  then 
* • T A 


STEP  3 
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If 


«3  > 1 + Cw 
If 

^3  < " ew 


If 


select  new  points  6 and  7 upstream 
one  point  and  restart  STEP  3. 

select  new  points  6 and  7 downstream 
one  point  and  restart  STEP  3. 


- ew  < fl3  £l  + ew  proceed. 

r3  - rfi  + q3(r?  - r6> 


STEP  4 


2(JX  - J2) 

l3  + l2  + g5  + g3  + Cg^  (q3  + q2  + h5  + h3) 


v ■ u C 
83  83  83 


P - J,  - 
1 


u (G.  + G,)  + v (H,  + H,) 
g3  5 3 g3  5 3 

~ - 

(S,  + Sc)  {?  - ) 


go  g 


+ 


3 y w g3  g5 


Steps  2 through  4 are  iterated  until  ve  calculated  variables  change 
by  less  than  a prescribed  amount,  or  until  a uia::imum  number  of  iterations 
is  reached. 
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5.6.47  Subroutine  XS1.P 

This  subroutine  Is  identical  to  subroutine  SLP  (Section  5.5.15). 


5*7  Link  60  Subroutines 


This  overlay  and  subroutines  contain  the  TBL  Module.  The  subroutine 
write-ups  are  from  Reference  5 with  few  exceptions. 

Only  a minimum  of  modifications  have  been  made  to  the  TBL  module.  These 
modifications  deal  mainly  with  the  linkages  between  the  various  modules  of  the 
SPP  code,  0DE  and  TD2P,  and  the  TBL  module. 


5.7.1  Program  Link  60 

This  subprogram  is  mainly  intended  to  facilitate  the  conversion  of  the 
SPP  code  overlay  structure  to  other  computer  systems.  The  only  function  of  this 
routine  is  to  call  subroutine  TBL. 
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5.7,2  Subroutine  BARCON 


Subroutine  BARCON  accepts  initial  conditions  and  executes  a Runge-Kutta 

solution  for  the  boundary  layer  along  the  contour.  The  initial  conditions  are  set: 

0=  0Q,  e = 90,  at  x=x0.  A step-size  Ax,  used  to  increment  the  axial  coordinate  x, 

is  determined  from  input  quantity  Axw=v  and  the  local  entires  in  the  x-table.  Having 

in  ax  i i 

two  first  order  ordinary  differential  equations  dg/dx  = f(x,  0,  0)  and  d0/dx  = g(x,9,<fl) 
and  using  subroutine  BARPRO  to  evaluate  the  derivatives,  a four  term  Runge-Kutta 
numerical  solution  is  used. 

1.  Evaluate 

fi  = _a£l  ' 9i  = !xl 

x^  x„ 

o o 

Set 

9 = 9o  + T"  (flJ  ' $ = + fep 

2.  Evaluate 

f2  = ^Xo+Ax.  'g2“dFlXo  + ^c 
Set 

Q=0o  + “^r^V  # + T”  k2) 


3.  Evaluate 


Set 


0-  + Ax  (f3)  , 0 = 0O  + Ax(g3) 


4.  Evaluate 


f 


_ d0  | 
4 "3x ' 


x. 


+ ax 


+ Ax 


5.  At  x + ax,  0and  0 are  then  evaluated  using 

0 = 9o+if  (f  1 + 2f2  + 2f3  + f4J 
0>  = *0  + (gl  + 2g2  + 2g3  = 94) 
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Now  using  these  values  of  x,  $ , 0 as  the  Initial  values  for  the  next  Interval,  we 
repeat  the  above  procedure  on  through  to  the  end  of  the  defined  contour. 

This  subroutine  also  calculates  and  prints  a table  of  normalized  contour 
points  corrected  for  displacement  thickness.  This  print  Is  made  to  assist  in  the 
preparation  of  a potential  flow  nozzle  contour  table  for  input  ot  the  TD2P,  TDK, 
or  other  method  of  characteristic  programs.  This  dimensionless  potential  flow 
(adjusted)  nozzle  contour  Is  calculated  by  the  following  relationships. 

POTENTIAL  = yt  * yt,REAL  " 5t 


^POTENTIAL  = [y  * yt(REAL  - 6*  cos  or]/y^ POTENTIAL  - 

Potential  = ' yt,  real  + 6*  sln  « 3/ytj  potential 

where 


yt,REAL  - actual  (geometrical)  throat  radius 
x are  XI  TAB 
y are  YITAB 


oi=  arc  tan  dy/dx. 

The  nozzle  throat  is  indicated  by  the  subscript  t. 
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5.7.3  Subroutine  BARPRO 


Subroutine  BARPRO  calculates  the  boundary  layer,  given  energy  thickness 
(0)  and  momentum  thickness  (0)  and  d0/dx  and  d0/dx  along  the  contour  at  a point 
(x).  The  interpolation  routine  XNTERP  is  used  to  evaluate  inviscid  flow  and  contour 
properties  and  derivatives  at  point  x.  Subroutine  ZETAIT  is  then  entered,  and  the 
parameter  f and  boundary  layer  thickness  A/  6#  6*  are  computed.  An  iteration  pro- 
cedure is  then  used  to  calculate  skin  friction  coefficient  (C^)  and  Stanton  number 

<°B>- 


1. 

2. 


An  initial  guess  (C,a  ) is  made 

f g 

T 

Calculate  CfR-  = (^L)  X"m  Cfa  R0 


3. 

4. 


5. 


Subroutine  CFEVAL  is  used  to  evaluate  Cf  as  a function  of  C^R^- 


Calculate  (-y— ) and  CfQ  = — 


'i 


__  T m 

aw  ^ aw  j / s j 


— j — ) 

e aw 


Cfa-Cfa 


A relative  error  comparison  Is  made.  If  | — ^ — — 2— | < TOLCFA 

fag 

convergence  is  considered  to  be  attained.  If  unconverged,  a form  of 
Wegstein's  method  is  used  to  calculate  a new  guess  (C^  ) and  steps  2 
through  5 are  repeated  up  to  a maximum  of  50  iterations.  9 


The  derivatives  d0/dx  and  d0/dx  are  now  evaluated,  from  the  differential 
equations,  at  the  point  x.  For  the  point  x at  the  end  of  a Runge-Kutta  step  (IND=1), 
total  heat  transfer  and  drag  quantities  are  also  computed. 


This  subroutine  also  prints  the  program  output  under  control  of  the  input 
variable,  IPRINT. 

This  subroutine  was  also  modified  to  print  and  write  (punch)  out  the 
linkage  data  needed  by  the  SPP  master  control  module.  In  particular  this  routine 
calculates  the  force  decrement  due  to  a turbulent  boundary  layer  using  the  equation 


AF 


(2iry)k  oe 


0 • COS  Cl  (1 


e 


) 


-1  dv 

where:  <*  = tan  -gj- 

K = 0 for  planar  flow 
K = 1 for  axlsymetric  flow 

and  the  decrement  in  I is  calculated  from 

sp 

^sp ' aF/l" 

where:  m = mass  flow  from  input  or  the  TD2P  module 

= force  decrement  as  calculated  above 
The  above  items  are  only  printed  out:  when  the  wall  slope  is  positive. 
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5,7,4  Subroutine  BARSET 


This  subroutine  coinputes  program  constants  and  sets  up  C and  H versus 

P 

T tables  and  tables  of  P and  T versus  x for  the  inviscld  flow.  Based  on  the  input 
table  of  C versus  T,  subroutine  BMTAB  is  used  to  generate  coefficients  of  a 

r 

cubic  spline  fit.  These  coefficients  are  then  used  to  compute  piecewise  partial 
integral  values  to  be  used  by  subroutine  SEVAL  which  relates  C^,  T,  H,  S and  P. 
Subroutine  GETPT  is  then  used  to  evaluate  pressure  and  temperature  at  each  point 
of  the  Mach  number  table. 


5.7.5  Subroutine  BMTAB 

This  subroutine  computes  the  coefficients  of  a piecewise  cubic  spline  fit 
of  tabular  data , 


5.7.6  Subroutine  BMFIT 

This  subroutine  computes  the  coefficients  of  a piecewise  cubic  spline  fit 
of  tabular  data. 
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5,7.7  Subroutine  CFEVAL 


Subroutine  CFEVAL  evaluates  the  adiabatic  skin  friction  coefficient,  C^, 
as  a function  of  C^R-,  based  on  Coles'  relationship.  For  values  of  C^R— , such 
that  0 <CfRg  <2.51,  Cf  is  evaluated  from 

~ .009896 

f (cfv*562 

For  values  of  C^R— , such  that  2.51  < C^R— < 1982759.2,  C^  is  evaluated  from  a 
piecewise  cubic  spline  fit  of  log  (C^)  vs  log(CfR^).  All  other  values  of  C^ 
are  considered  out  of  range. 


5.7.8  Subroutine  DIRECT 

Subroutine  DIRECT  controls  the  overall  flow  of  the  program.  It  calls  sub- 
routines which  successively  read  in  data,  compute  program  constants,  fit  curves 
and  then  solve  the  boundary  layer  solution  algorithm. 


5.7.9  Subroutine  EDITCP 

This  subroutine  attempts  to  edit  a table  of  Cp  versus  T from  the  boundary  layer 
edge  condition  table  of  Cp  versus  T if  the  IFEDGE=1  or  2 option  was  selected  in 
the  $TBL  namelist  and  if  n<?  C versus  T table  was  input. 

r 
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5f7 . 10  Subroutine  FIIF 

Subroutine  FIFF  defines  the  function  (of  "s")  to  evaluated  by  the  numerical 
integration  routine  INTZET.  Three  second  degree  coefficients  (Ap._Bp,  Cp)  and  an 
exponent  value  (mf)  are  set  by  subroutine  ZETAIT.  Enthalpy  value  h is  computed 

from 

h = Ap  + Bps  + CpS® 

Entering  subroutine  SEVAL  with  h.  we  obtain!.  The  appropriate  form  of  the  function 
is  then  evaluated: 

For  IP  = 1,  f(s)  = (-r)  s mf  (1  - s) 

F t 

I„  - 2.  f(s)  - (il  s »« 
f t 

where  T/r  has  been  used  for  q/q  . 
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5.7.11  Subroutine  GETPT 


For  a given  Mach  number,  subroutine  GETPT  computes  pressure  P and 
temperature  Te  from  . H and  T relations.  For  constant  specific  heat  (y  = yQ ) 


i y “ 1 

(1 


For  variable  Cp  the  following  iteration  is  used 

1.  An  initial  guess  is  made  for  temperature  T from 

g 


T 

g 


M2 


2.  Using  subroutine  SEVAL,  Cp  and  H are  computed  at  temperature  T 

3.  The  specific  heat  ratio  y is  then  evaluated  from 


C 


4.  Temperature  T is  then  calculated 

2(H  - H) 

o 

e y RM2 

5.  A relative  error  comparison  is  then  made.  If 


T - T 

!-^P — 2-'  < TOLZME 

g 

convergence  is  considered  to  be  attainea  and  subroutine  SEVAL  is  used  to  evaluate 
pressure  P0  at  temperature  Tg.  If  the  convergence  criterion  is  not  satisfied,  a 
form  of  Wegstein's  method  is  used  to  calculate  a new  guess  for  Tg  and  steps  2 
through  5 are  repeated  up  to  a maximum  count  of  50  iterations. 
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5.7.12  Subroutine  INTZET 


Subroutine  INTZET,  given  a function  definition  and  limits  of  Integration 
(x^  to  Xg),  uses  a numerical  technique  based  on  the  qulntlc  spline  to  evaluate  the 
Integral  of  the  function  from  to  x2 . The  Interval  over  which  the  Integration  is 
to  be  performed  Is  divided  Into  20  coarse  subintervals.  Since  the  limits  are  re- 
stricted to  certain  values  (0,  1,  £ , l/£ ),  a check  on  the  length  of  Interval  (x2  - 
Xj)  Is  made.  If  (x2  - Xj)  < .10,  the  function  Is  evaluated  at  the  21  sublnter/al 
endpoints  and  the  Integral  Is  approximated  by  the  qulntlc  technique,  with  a "one- 
sided Simpson's  Rule"  being  used  over  the  extreme  subintervals. 

For  (x2  - Xj)  > .10,  the  maximum  value  of  the  function  within  the  Interval 
Is  found,  and  the  approximation  is  performed  over  the  coarse  steps  up  to  a point 
where  the  value  of  the  function  Is  20%  of  the  maximum  value.  Beyond  this  point, 
the  remaining  coarse  sublntexvals  are  each  further  divided  Into  3 fine  subintervals. 
The  function  Is  evaluated  at  the  endpoints  of  these  fine  subintervals.  The  Integral 
approximation  Is  completed  using  this  finer  step. 

For  a given  function  f(x)  evaluated  at  equally  spaced  (ax)  points  xQ,  x^, 
x2'  x3#  the  qulntlc  approximation  for  Integration  Is  given  by: 

J 2 fW  d.  = (-f(x0)  + 13f(xj)  + 

13f (x2)  - f(x3)) 

The  "one-sided  Simpson's  Rule"  is  given  by: 
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5,7,13  Subroutine  QUITS 


Subroutine  QUITS  is  entered  when  an  error  has  been  detected  in  the  input 
or  if  an  unreasonable  quantity  has  been  detected  during  execution.  It  prints  out 
appropriate  error  statements  and  the  contents  of  labeled  common  storage. 


5.7,14  Subroutine  READSI 

Subroutine  READSI  is  a modified  version  of  subroutine  READIN  of  Reference 
5.  The  name  was  changed  to  avoid  a naming  conflict  with  a subroutine  by  the 
same  name  in  the  ballistics  module. 

This  routine  reads  the  input  data  for  the  TBL  module  in  the  $TBL  namelist 
as  described  in  Volume  III,  Section  2.10o  Subroutine  READSI  also  checks  the  in- 
put data  tor  consistancy  and  prints  out  dianogstics. 

In  addition  to  the  above,  this  routine  traps  and/or  reads  any  linkage  data 
generated  by  the  0DE  or  TD2P  modules. 
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5.7.15  Subroutine  SEVAL 

Subroutine  SEVAL,  using  certain  defined  Integral  relationships,  and  a piece- 
wise  cubic  curve  fit  of  specific  heat  (Cp)  as  a function  of  temperature  (T),  for  a 
given  Indicator  value  (Ind),  evaluates  the  appropriate  quantity:  Pressure  (P),  en- 
tropy (S);  enthalpy  (H);  temperature  (T);  specific  heat  (Cp) . 

For  Ind  = -1:  P = P(T,S) 

0:  S = S(T,P) 

1:  H = H(T),  Cp  = Cp(t) 

2:  T = T(H),  Cp  = Cp(t) 

3:  T = T(S,P) 

The  basic  relations  are: 

Cp  = A + BT  + CP  + DT3 
H(T)  = H.  + fT  C (t)  dt 

1 md  P 

T CnW  P 

S(T)  = S(Tq)  + JT  — E — dt  - R In  (y-) 

T o 


5-209 


1 


5. 7 ,16  Subroutine  START 

Subroutine  START  locates  the  sonic  point  (M=l)  from  the  input  Mach  num- 
ber and  x tables  and  then  computes  initial  values  for  momentum  and  energy  thick- 
nesses. 


The  table  of  Mach  numbers  is  first  scanned  to  locate  an  exact  value  of 
M=l.  If  none  is  found,  the  sonic  point  location  is  determined  from  the  following 
iteration  procedure. 

1.  Locate  the  i-th  interval,  in  the  table,  which  brackets  M=1  and  initially 
guess  xg=xi# 

2.  Using  subroutine  XNTERP,  evaluate  Mg  at  xg. 

3.  An  absolute  error  comparison  is  made.  If 

[ Mg  - l.|  < TOLZME 

convergence  is  considered  to  be  attained  and  appropriate  flow  quantities 
are  evaluated  at  x . If  unconverged,  a secant  method  is  used  to  cal- 
culate a new  guess  for  x and  steps  2 and  3 are  repeated  up  to  a max- 
imum of  50  iterations.  g 

Using  the  flow  quantities  just  computed,  subroutine  INTZET  is  used  to 
evaluate  integrals  I ^ and  I2  for  £=  1.  The  shape  factor  5*/0  is  then  evaluated 
from  equation  (125  of  Ref,  5).  An  iteration  for  the  skin  friction  coefficient  Cf 
is  entered  (similar  to  that  in  subroutine  BARPRO),  and  a value  of  momentum  thick- 
ness q is  computed.  This  value  is  assumed  to  be  the  initial  value  for  0and  </, 

at  x, . 

3 
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5.7.17  Subroutine  TBL 

Subroutine  TBL  initiates  the  execution  of  the  TBL  module  by  calling  sub- 
routine DIRECT. 


5.7.18  Subroutine  XNTERP 

Subroutine  XNTERP,  given  tables  of  n points  of  independent  variable  (x^) 
and  dependent  variable  (y^,  returns  point  value  (y)  and  first  derivative  (y1)  for  a 
given  value  of  x. 


A local  quintic  spline  interpolation  is  defined  as  follows:  The  independent 
variable  table  is  searched  until  the  points  bracketing  x are  found  (x^  < x < *1+1). 
Using  the  four  surrounding  points  (x^,  x^  x[+1,  xl+2)#  coefficients  of  the 
parabola  which  pastes  through  the  three  leftmost  points  and  that  which  passes 
through  the  three  rightmost  points  are  obtained.  For  the  given  value  x, 

yL  = c1  + c2  (x  - x^j)  + c3  (x  - x^)2 

yR  = c4  +c5  (x  - x^)  +c6  (x  - xf 

?L  = c2  + 2c3  (x  - 

?R  = c5  + 2c6  Oc  - xt) 

A cubic  weighting  function  ^and  its  derivative  a'  are  defined  as  follows: 


U(x) 


x - x. 
xi+l  “ xi 


ot  (x)  = 3U2  - 2 IP 
' (x)=  (6U  - 6U2)  U'  (x) 


Vhe  dependent  variable  y and  its  derivative  y'  can  now  be  evaluated. 
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5.7.19  Subroutine  ZETAIT 


Subroutine  ZETAIT  calculates  the  shape  parameter  £ and  boundary  layer 
thicknesses  a,  fi  and  5*  at  a point  x,  for  given  values  of  0 and  0 . 


For  £ > 1 


£ = (i  //n~I2'“13 
0 1 I, 


■) 


where 


For  £<1 


6 = 


C = 


6 


± 1 
e i^q7ic 


i 

n+l' 


* l/n  I6  I? 

T +T*1"  } 

0 a4  a5 


6 = 


where 


C = 


± !& 

e i- 


n+l 


An  iteration  procedure  is  used  to  evaluate  £. 

1.  An  initial  guess  r is  made 

9 

2.  For  the  value  of  £ ^1  or  < 1)  appropriate  functions  of  p/p  have  been 

defined,  and  thr  proper  integrals  are  evaluated  using  subroutine  INTZET 

3.  £ is  calculated  by  the  appropriate  formula  a&ove 

4.  A relative  error  comparison  is  made.  If 


£ - £ 

| — t—2-I  < TOLZET 

C g 

convergence  is  considered  to  be  attained,  if  not  converged,  a form 
of  Wegstein's  method  is  useo  to  calculate  a new  guess  for  £ and 
steps  2 to  4 are  repeated  up  to  a maximum  of  50  Iterations.  g 

The  boundary  layer  thicknesses  5*  and  A are  then  calculated. 
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Appendix  A - Conversions  Aides 


The  SPP  code  was  developed  under  the  rules  and  limitations  of  the  SCOPE 
3.3  operating  systems  of  CDC  6000  computers.  The  following  notes  are  intended 
to  aid  the  conversion  of  the  SPP  code  to  other  computers  and/or  operating  systems. 

1)  Subroutine  FIND  of  Section  5.1.6  is  not  the  same  routine  as  described 
in  Section  5.6.5.  On  operating  systems  which  poses  a problem. 

It  Is  suggested  that  the  name  of  the  routine  of  5.6.5  be  changed  to 

FINE.  This  routine  is  referred  to  in  subroutines  AC0MP,  0NED,  PARTTL, 
and  TRACE. 

2)  Subroutine  ERROR  of  Section  5.4.8  is  not  the  same  as  the  routine  de- 
scribed in  5.6.37.  On  operating  systems  where  this  poses  a problem 
It  is  suggested  that  the  routine  of  5.4.8  should  be  changed  to  the  name 
(external  reference)  of  ERR0S.  This  routine  (5.4.8)  is  only  referred  to 
in  subroutine  RE  ADI  N. 

3)  A possible  conflict  in  labeled  common  areas  occurs  in  LINK41  and  LINK50 
routines.  Labeled  common  blocks  by  the  name  NAMEW  occur  in  both 

of  these  overlays,  however,  these  are  not  intended  to  be  the  same.  It 
is  suggested  that  the  common  block  in  subroutine  0DKENP  be  changed 
to  the  name  NAMW  to  avoid  this  possible  conflict. 

4)  Contrary  to  the  rules  of  F0RTRAN  IV,  EQUIVALENCE  statements  can 
sometime  increase  the  size  of  labeled  common  areas.  This  can  cause 
problems  on  some  operating  systems.  In  particular  the  labeled  common 
block  FNALBT  in  subroutine  DERIV  can  be  made  to  be  4 locations  (cells) 
longer  than  defined  in  subroutine  MAIN1D.  If  this  causes  problems, 
the  length  of  this  common  block  should  be  extended  by  4 dummy  loca- 
tions (which  are  not  used)  in  subroutines  MAIN  ID,  FLU,  and  IAUX. 

Other  minor  corrections  to  the  SPP  code  which  will  update  the  SPP  code  to 
the  1.2  version  are: 

a)  Moving  the  labeled  common  block  D0L00PS  into  the  0DK  routine. 

b)  Adding  the  variable  CPS0L  to  the  $TD2  namelist. 

c)  Adding  the  test  at  approximately  line  65  of  subroutine  TD2P  of 

IF  (FTD2P  .EQ.  0.0)  CPS0L  = 0.0 

d)  Changing  the  test  is  subroutine  AGP  from 

IF  (FTD2P  .EQ.0)  G0  T0  12 
to:  IF(CPS0L  .EQ.0)  G0  T0  12 
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The  conversion  of  the  SPP  code  to  the  Univac  1108,  EXEC  8 operation  system 
took  approximately  one  man  day  of  effort.  However,  it  is  recognized  by  long  and 
painful  years  of  experience  that  the  authors  of  a particular  code  have  much  less 
trouble  converting  codes  than  those  unfamiliar  with  the  particular  program  in 
question.  Therefore,  questions  concerning  the  SPP  code  should  be  directed  to  the 
Air  Force  contracting  officier  in  charge  of  the  SPP  code  and  responses  will  be  made 
through  him. 

Other  notes  concerning  conversions  of  the  SPP  code  are  as  follows  for 
Univac  and  IBM  users. 

I)  The  cards  with  0VERLAY  (SPP,i,j)  should  be  removed. 

II)  The  cards  with  PR0GRAM  LINKij  should  be  changed  to  SUBR0UHNE 
LINKiJ  and  a RETURN  card  should  be  added  just  before  the  END  card. 

III)  The  calls  to  0VERLAY  (3LSPP,l,j)  should  be  changed  to  CALL  LINKij. 

It  is  hoped  that  these  brief  notes  will  substantially  aid  in  the  conversion  of 
the  SPP  code  to  other  computer  and/or  operating  systems  for  which  the  SPP  code 
was  not  designed  for. 
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